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© A technique for the detection and removal of local defects in digital continuous-tone (mages. 
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© The present invention is a method for automatically detecting and correcting a wide range of local digital 
image defects with minimal user intervention. The detection process employs brightness and color thresholds in 
conjunction with magnitude thresholds on residuals of nonlinear spatial filters to separate defects from scene 
content with minimal confusion. The detected defects are then cosmetically corrected by combinations of 
nonlinear smoothing and grey-scale erosion. Several options are outlined for the feature selection, detection, and 
cleaning operations depending on source type and computational constraints. 
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Cross-reference to Related Applications 

The present application is related to U. S. Patent Application Serial No. 07/934,089 entitled "Process for 
Detecting and Mapping Dirt on The Surface of a Photographic Element" filed 08/21/92 by Robert Gray et 
5 al., corres. to EP 0 585 759 A1. 

Field of the Invention 

The present invention is related to the field of scanning images, primarily from photographic film, to 
10 form electrical equivalents of the scanned images and more particularly to the detection and the removal of 
defects due, for example, to dirt contamination or physical damage of the film. 

A portion of the disclosure of this patent document contains material which is subject to copyright 
protection. The copyright owner has no objection to the facsimile reproduction by any one of the patent 
disclosure, as it appears in the Patent and Trademark Office patent files or records, but otherwise reserves 
75 all copyright rights whatsoever. 

Background of the Invention 

Digital images created from electronic scanning of continuous-tone photographic film often reveal 

20 visually objectionable defects due to dirt contamination or surface damage of the film. These defects 
usually must be digitally corrected by manual digital retouching (e.g. via "dustbu sting" and other cloning 
techniques), which require that the operator visually identify each local defect in the image. 

A patent of interest for its teaching in this art is U.S. Patent No. 4,189,235. entitled "Test Device For 
Dynamically Measuring The Degree of Dirt Accumulation On Bank-Notes" by Guter et at. That patent 

25 describes a method for inspecting opaque web materials for dirt accumulation. A reflected signal from the 
web material is sensed by three adjacent photosensors as the material moves past the sensors. When the 
signal sensed by the center sensor is significantly different than from that sensed by the side sensors the 
sensed area is classified as containing dirt. 

In conventional optical printing systems, if dirt or a scratch appears on the photographic image, the 

30 photographic negative or slide may be cleaned and the print remade or the print itself may be retouched. In 
the case where the images are written to an optical disc, it is not convenient nor desirable to rewrite the 
image, because many images are written onto a disc before the images are viewed, and once an image is 
written it cannot be removed. Therefore, it is desirable to monitor the dirt present in photofinishing 
environments and to assess the effectiveness of any selected film cleaning methods. 

35 tn U.S. Patent No. 4,907,156, entitled "Method and System for Enhancement and Detection of Abnormal 
Anatomic Regions in a Digital Image" by K. Doi et al. there is disclosed a technique for taking the 
difference between an "enhanced" and a "degraded" version of the input image to identify and remove 
structured anatomic background. A series of single thresholds are then applied to the difference images 
using the shape/size behavior of image regions to perform a feature extraction that is related to the 

40 detection of lung nodules. The feature extraction is implemented by the use of thresholds against which 
difference image pixels are compared to identify which super threshold "blobs" are lung nodules and which 
are not. 

A significant problem in the detection of dust is to provide an adequate way of detecting "small" 
objects. There are several possible approaches to this problem. For example, one could color-classify all 

45 the color pixels in an image, then spatially connect the pixels of a common color-class into segments, and 
then examine the spatial dimensions of the segment. This approach, while possible, has the drawback that 
the color classification must be quite accurate, or else some marginally-colored pixels will be misclassified. 
The problem is that differences in the direct RGB (or luminance-chrominance) image values do not alone 
give much contrast information between, say, a white small dust artifact and a bright grey extended scene 

so region on which the artifact is superimposed. Thus after color-classification one may be left with very large 
segments that are reasonably neutral and bright, but which now mask visible dust artifacts. 

A better approach prior to classification is to preselect only those points which have large local contrast 
from their surround. This may be directly accomplished by (1) spatially smoothing the color image and (2) 
subtracting the smoothed color image from the original color image. The resulting difference (or "residuals") 

55 image will emphasize all color pixels which have a large magnitude compared to some weighted local 
average magnitude. The key decision here is then the choice of spatial smoother. The most obvious first 
candidate is a low-pass linear filter (e.g. a simple spatially weighted averaging of all grey levels (GLs) in a 
local window); the properties may be easily analyzed mathematically, and (perhaps more important) such 
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fitters may be rapidly implemented in software and hardware (the residual image is basically a high-pass- 
filtered version of the original). Unfortunately, linear residuals images have the property of retaining all local 
edges in an image, including very long edges due to scene discontinuities. The result is a very "cluttered" 
image (i.e. the residuals image contains many bright pixels having no connection with dirt or scratches); a 

5 significant amount of additional time must then be spent in scene reasoning to clean up the residuals map 
before defect decisions can be made. In short, linear-filter residuals poorly discriminate between local 
defects and true scene details. 

The identification of anomalies such as dust or scratches leads to the desirability of removing the same 
without the intervention of an operator. The present invention is directed to a methodology for identifying 

10 and cosmetically correcting such anomalies with minimal operator intervention. 

Summary of The Invention 

One of the preferred methods of the present invention is: A method for the detection and removal of 
75 defects in digital images comprising the steps of: 

a) preprocessing a digital image represented by pixel values to form a feature image in which defect 
pixels have large local contrast from their neighboring non-defect pixels; 

b) testing the value of each feature image pixel value to determine if the value is within a range of pixel 
values expected, of defects; 

20 c) edge-preserving spatial filtering of each of the feature image pixel values that are within the range of 
expected defects; 

d) forming a residual value for each pixel as a function of the difference between the corresponding 
feature image pixel value and the filtered value of step c); 

e) testing each residual value of step d) to determine if the residual value is within a range of residuals 
25 values expected of defects; 

f) forming a map in which each pixel that is within the range of expected defects according to steps b) 
and e) is marked as a defect; and 

g) correcting the digital image using the map created in step f). 

The present method, is called an autodustbuster (ADB) algorithm, and it has been successfully tested 
30 on both consumer-type color photographic images and color motion picture frames with good results. 

From the above it can be seen that it is a primary object of the present invention to provide an 
automatic local defect detection and removal technique for digital images. 

The above and other objects of the present invention will become more apparent when taken in 
conjunction with the following description and drawings wherein like characters indicate like parts and which 
35 drawings form a part of the present invention. 

Brief Description of The Drawings 

Figure 1 illustrates a system on which the present invention may be implemented. 
40 Figure 2 illustrates the windowing of an image raster to determine the grey level of a central pixel GLc 
at a position X. 

Figure 3 illustrates in chart form a histogram of GL values for an associated window. 
Figure 4 illustrates the averaging of the samples "nearest" in grey level GL to the central pixel GLc for 
the histogram of Figure 3. 

45 Figure 5 illustrates the replacement of the GLc value with the grey level average GL av at the position X. 
Figure 6 illustrates the process flow for a multi-band original image. 
Figure 7 illustrates the process flow for a resolution pyramid of an image. 
Figure 8 illustrates the process flow as applied to an R, G, B image. 
Figure 9 illustrates a flow for ADB processing within a resolution pyramid. 

50 

Detailed Description of the Invention 

Referring to Figure 1, the present invention is implemented using a digitizing scanner 10, a control and 
logic device 24, such as a computer with a display device 26, and an optical disc writer 36. As shown, a 
55 strip of film 15, containing one or more frames 18 of developed film, is placed into a scan gate 12 where 
the strip of film 15 is scanned under direction of control and logic circuitry 30. As each frame 18 is scanned 
the resultant scanned image data, represented by block 22, is digitized and transmitted to a memory 28 for 
storage. The computer 24 process the stored data, in a manner to be described, to provide output image 
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data 32 which may be written to an optical disc 34 by the optical disc writer 36 to provide a report as to the 
characteristics of the anomalies. The scanning device 10 is capable of quantizing pixel values into multiple 
brightness levels in separate red, green, and blue channels. A minimum number of brightness levels would 
be approximately 64 for an adequate quality image with a typical number being 256. Although the image 

5 data 22 has been described as originating from the operation of scanning film, it is wefl understood that 
other well-known techniques may be used to provide the image data. 

The autodustbuster (ADB) algorithm and several processing options that allow customization for 
different source types and computational constraints will now be described. 

The required ADB operations may be divided into (1) feature selection, (2) defect detection, and (3) 

w defect correction. The operations must ensure that the maximum number of defects is detected and 
optimally corrected while minimally affecting scene elements which are not defects. Thus, the more 
confidence which the user has in the accuracy of the detection process, the more extreme can be the 
removal process; conversely, as one's confidence in avoiding "false alarms" from scene elements 
diminishes, milder removal methods may be used to avoid damaging real scene elements. 

75 Visible dirt artifacts in digital positives, created from scanned color negatives, possesses the following 
properties: they are "white" (low color saturation), "bright" (very high apparent exposure values), "small" 
(small spatial extent) and "contrast" (large differences from their immediate surround). (If dirty color 
positives were instead scanned, then the dirt artifacts would be "dark" rather than "bright"). In addition, 
many such fragments exhibit little brightness variation within individual defects. 

20 It is obvious that the above characteristics are rather loosely defined, and that many real scene 
characteristics may possess one or more of these attributes, Thus the goals of an ADB algorithm are to 
maximize the percentage of defects detected ("hits") while minimizing the number of real scene elements 
mistakenly considered to be defects ("false alarms"). These goals tend to be competitive, that is, setting 
thresholds on features to increase the percentage of "hits" tends to increase the number of "false alarms" 

25 and vice versa. A detection mechanism must be used in which defects appear as dissimilar as possible 
from true scene elements. 

In most imagery, spatially small objects which differ greatly in brightness or color from the immediate 
surround are relatively rare. The ADB algorithm exploits this property by the innovative use of a particular 
form of nonlinear spatial filter. The properties of this filter are described in the next section. 

30 Figure 2, illustrates in flow chart form the generic autodustbuster process flow. The original image of 
block 100 is a multi-band digital image represented by one or more matrices of pixel values or equivalents 
grey levels (GL). In preprocessing block 1 1 0 some effective and convenient remapping of the original image 
may be performed. An example of preprocessing is the digital resampling of the image to reduce its size 
for computational expediency. An alternate preprocessing is the selection or recombination of the original 

35 image bands to produce a feature image with fewer bands. The detection suitability of the resulting feature 
image, in any case, should be at least comparable to that of the original image. Oval 120 represents the 
feature image that is outputted from block 110. The feature image is subjected to a band brightness test in 
block 130. The test applies low and high threshold evaluations to every pixel of each band in the feature 
image. Pixel values that are above the high threshold or below the low threshold in each band are 

40 considered to have passed the test and are forwarded to an edge-preserving-smoothing (EPS) filter 140, 
This filter non-Iinearly smoothes the selected pixel in all bands in a manner that retains major scene edges 
as subsequently described. The output of the filter at this pixel position is written to a corresponding 
position in the smoothed image 142. For those pixel values that are below the high threshold value, but 
above the low threshold value no subsequent EPS filtering is applied, i.e. the test is considered failed. Each 

45 band of the feature image pixel is copied to the resultant output of the smoothed image 142 and does not 
undergo EPS filtering. At difference node 145 the pixel values of each band of the smoothed image are 
subtracted from the pixel values at the corresponding position in the feature image 120 resulting in the 
residuals image 148. 

All of the pixel values from the residuals image are subjected to the residuals test of block 150. The 
so residuals test consists of a minimum and a maximum threshold that is applied to each band of the residuals 
image 148. When bright defects are to be detected, for a pixel to pass all band values must exceed the 
minimum threshold and at least one band value must exceed the maximum threshold. When a dark defect 
is to be detected a pixel will pass if all band values are less than the maximum threshold and at least one 
band value is less than the minimum threshold. 
55 Pixels that pass are forwarded to a classmap 160 and marked with a distinctive value. The classmap 
160 is a single band image representation which is used to mark the location of detected defects. Pixels 
that fail are marked with a different value in the classmap 160. In a classmap conditioning block 170 the 
values within the classmap 160 may be further modified to accommodate known spatial characteristics of 
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the defects. For example, a one-pixel morphological dilation of the defect-marked (passed) pixels in the 
classmap is typically performed. The conditioned classmap 175 is directed to the cleaning block 180 long 
with the original image pixel values. Pixels in the original image that are not marked as passed in the 
classmap will be copied unchanged to the output cleaned image 190. Pixels which are marked as passed, 
i.e., defect, are corrected using neighboring non-defect pixel values. The corrected values are then written 
to the cleaned image 190. 

A more detailed discussion of the major blocks of Figure 2 will now be undertaken. In ADB processing 
the feature image 120 used to detect the defect locations may be different from the original image 100 
which will be corrected. This is generally advantageous when computational time is at a premium and EPS 
filtering 140 of the full (generally three-color) image is too burdensome. In this case it is often possible to 
apply a preprocessing module 110 to create a one- or two-band "feature" image 120 which instead will be 
EPS-filtered and classified for defects. The effectiveness of this strategy will depend on the success with 
which a reduced-dimensionality feature image will retain the size, brightness, contrast, and color defect 
discrimination of the original image 100. This will of course be application-specific. For example, in digitized 
imagery originating from superimposed color cells (as in animated motion pictures) defects may be of any 
color and saturation; thus little is lost in detection performance using a reduced-dimension space which 
sacrifices color and saturation information. In this case good success has been achieved in using a "band- 
, maximum" one-band feature image 120 band where every feature pixel value is the maximum value of the 
input red, green, and blue bands. By contrast, in Photo CD applications the principle source of defects is 
dirt on the negative during scanning, which appears as "white" specks in the digital positive; here saturation 
information should be retained as an important identifier, and two- or three-band EPS-filtered feature images 
are the norm. Figure 7 illustrates the retention of all three color original color bands for classification. 

Rescaling of the feature image grey-level range during preprocessing 110 is often advantageous. In 
current implementations this consists of linearly scaling to a 0-255 pixel value range even for 10- or 12-bit 
original data. This results in a slight speed improvement during EPS-filtering 140 when the histogram is 
searched.at every window location (see Figures 4 and 5 along with associated discussions). In Photo CD 
applications (e.g. Figure 7 this rescaling may be determined by finding .5% and 99.5% cumulative 
histogram points of the color-corrected bands of the entire input image. The minimum .5%-point and 
maximum 99.5%-point from among all bands are mapped to 0 and 255 pixel values, respectively in the 
feature image, thus ensuring that no color shift occurs. If the input data is not color balanced it may be 
necessary to perform a color balancing step on the original data. 

EPS Filters and k-Nearest Neighbor (kMM) 

To detect the critical size of the desired defects a smoothing filter is desired which retains extended 
image edges but smoothes small local contrasty image regions (blobs). Such fitters exist as various 
nonlinear noise-reduction methods under the category of "edge-preserving smoothers" (EPS). The most 
widely known EPS filter is the median filter, which is one member of the class of rank-order filters. In rank- 
order filters, the pixel magnitudes are sorted and an average of some contiguous number of the sorted 
magnitudes is computed. For example, in a median filter the "middle" element of the sorted magnitudes is 
retained (an average over one sample). The sample mean is, at the other extreme, also a rank-order filter 
(although a poor EPS); it provides the average of all sorted magnitude values. 

A poor property of standard rank-order filters for image operations is that they are insensitive to the 
relative spatial location of the pixels within the current window; e.g. the magnitude chosen as the "median" 
may exist at the center of the window or at an extreme edge of the window. A number of specialized 
nonlinear window operators have-been devised to account for local spatial structure. One example is the 
sigma filter (J.-S. Lee, CVGIP,24,255-269. 1983). In this filter a histogram is created of the grey levels in a 
sliding window; based upon an estimate of the image noise, a range of histogram values about the grey 
level of the pixel at the center of the window is computed. One characteristic of the sigma filter, however, is 
that it is specifically designed for noise removal and its averaging range is solely dependent on the noise 
characteristics; the filter averages over a different number of pixels at each window location. For "blob" 
detection, however, what is desired is a filter whose parameters are sensitive to local shape characteristics. 

A relative of the sigma filter which gives this control is the k-nearest-neighbor (kNN) filter (Davis and 
Rosenfeld, IEEE Trans, SMC, 7, 107-109 ,1978). This filter is implemented in the EPS filter module 140 of 
Figure 2. This filter has the following operation per color band: 

Referring to Figure 3, slide a rectangular window of odd dimensions n x m over an image in raster 
order. At each location determine the grey level of the center pixel (GL C ) or equivalently pixel value, and 
also create a histogram of the grey levels in the local window as per Figure 4. Beginning at the histogram 
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bin containing the center grey level GL C , sum contents of the histogram bins on both sides of the starting 
bin in the following fashion: Test if the. contents in the bin for Gl_c are greater than or equal to k pixels. If 
this is true, then write out as the averaged value GL av just the existing value GL C . If the contents for bin GL C 
are less than k, then add the contents of the bin immediately below the bin and above GL C to the sum and, 
5 again test whether the sum is greater than or equal to k. If the sum equals k, then calculate the average 
grey level according to the formula 

GL av = SUM {GLi x hist(i)} 
w 

SUM {hist (i) } 

75 where the index i is over the summed histogram bins. If the sum is greater than k, then the extreme bins in 
the sum are proportionally reduced in weight in the calculation for GL av so that the total bin weight in the 
sum approximately equals k {as shown in the attached source code). If the sum is less than k, then the next 
bin below Gl_c and above GL C are summmed and tested against k. This process is repeated until the sum of 
histogram bins about GL C are greater than or equal to k. Continue until k pixels (k < n x m) have been 

20 summed. Find the average grey level of the summed pixels per Figure 5 and write this value to an output 
image at the location of the window center per Figure 6. 

The k-nearest neighbor filter has some interesting filtering properties which are not immediately obvious 
from the above description. First, any constant-GL image structure (blob) which has an area of k or more 
pixels in the local window will not be significantly smoothed; any structure of less than k pixels will be 

25 increasingly smoothed as k decreases. Second, at "ramp" edges between otherwise constant grey level 
regions the slope of the ramp will be increased (sharpened) by the filter. 

The first property is the one of significance for "blob" detection. For example, if the window size is 5 x 
5 pixels, and if k is, say, 15, then extended straight edges having a depth (constant pixel value in both 
directions perpendicular to the edge) greater than three pixels will be effectively unchanged, and thus 

30 disappear in a residuals image. "Blobs" of area less than 15 pixels that fit entirely in the window will be 
smoothed and thus will appear strongly in the residuals image. This is the desired residuals property for a 
local defect detector. 

A major concern for a k-nearest neighbor filtering is computational speed. This is not a linear filter and it 
is not separable, nor does it have other usable symmetry properties. The major operations per window 

35 location are histogram creation and then accumulation with testing of the total accumulated samples. The 
current k-nearest neighbor implementation does significantly reduce the cost of computing the histogram by 
merely updating the window histogram in a horizontal strip of image rather than performing complete 
calculations per window. This results in a fractional saving of approximately (n-2)/n in histogram-computing 
operations (for a n pixel x m line window), i.e. the savings increase as window length increases. (This fast 

40 histogram updating is equivalent to that used in the "fast-median" filter-Huang and Tang, IEEE Trans, 
ASSP, 27, 13-18, 1979). 

Residuals Test 

45 As mentioned above, a residuals image 148 is created as shown in Figure 2. The anomaly classifier 
steps contained in modules 150 and 170 accept as input the residuals image 148 and create as output the 
final classmap 175 which marks the pixels to be cleaned. 

The residuals test in module 150 creates an initial classmap by performing tests on each pixel of the 
residuals image 148. In the classmap conditioning module 170 a spatial segmentation (i.e. connected- 

50 component labeling) of the non-background pixels in the classmap may be optionally performed with 
additional tests applied to the same segment. 

In the simplest (single-state) version and for detecting bright defects, each color residual pixel is 
evaluated against two residual-GL thresholds, GL rmin and GL rmax . In this implementation, all color residual 
pixels whose minimum band residual GL exceeds GL rmin andwhose maximum band residual GL exceeds 

55 GL rmax are uniquely marked for subsequent cleaning; pixels which do not satisfy both of these constraints 
are marked to be left unaltered in the output image. (For single-band feature images GL rmin equals GL rmax ). 
The two-dimensional map of category markings per pixel is termed the "class map". These thresholds 
impose a minimum required contrast between a defect pixel and its immediate background. The (smaller) 
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GL rmin threshold is generally set slightly above the "zero point" in the residuals image (e.g. to 133 for 
rescaled features with residuals biased by a value of 128 as described above) and acknowledges that valid 
bright defect pixels , e.g. due to dirt on a scanned negative, should not have a negative contrast (residual 
values less than the zero point) in any band. The (larger) GL rma x threshold states a requirement of 
5 significant local contrast in at least one band for the defect to be visible. (Recalling that these thresholds are 
on the residuals GLs t it is obvious that not all bands should be required to pass a large GL rmax value; for 
example, a white dirt artifact against a bright red background would have large residual values in the green 
and blue bands, but only a small residual value in the red band.) 

For two-state anomaly classification, two sets of residuals thresholds are used, i.e. GL rmin i, GUmaxi. 
io GL rmin2 , and GL rmax2 . In general, the minimum-residuals thresholds are set equally, and slightly above the 
zero point as mentioned above; the maximum-residuals threshold is set somewhat higher for labeling high- 
confidence pixels than is the maximum-residuals threshold for the low-confidence pixels. This labeling 
results in a class map where the pixels possess one of three labels (normal, low-confidonce dirt, and high- 
confidence dirt). In two-state classification, this class map undergoes image segmentation in module 170, 
75 i.e. all residuals pixels spatially contiguous in the class map which have the same class are, assigned a 
unique code value. (These code values can be thought of as uniquely labeling all the "blobs" of suspected 
defects in the image.) A fairly efficient algorithm for segmenting the class map is employed which involves 
finding the lowest prototype code value in a spatial "tree" of prototype segments of the same class, and 
then renaming all the segments in the tree with this lowest code value. 
20 Following segmentation, all segments which contain low-confidence defect pixels and which are 
spatially touching segments containing high-confidence pixels are upgraded to high-confidence (i.e. the 
class map values of the pixels are reassigned to the "high confidence of an defect" class). The purpose of 
this step is to make use of the adjacent presence of very obvious defects to decide that a slightly less 
contrasty region is part of the same defect. The segments which remain at low confidence are either 
25 relabeled as "clean" (i.e. all pixels of the segment are marked for no cleaning in the classmap) or else 
name their pixels marked for subsequent cleaning state in module 180. The segments, which are marked as 
high confidence of being defects, have their pixels marked in the classmap for subsequent erosion cleaning 
in module 180. 

30 Cleaning 

Two general cleaning methods are presently available for use within module 180 in the ADB algorithms: 
EPS substitution and greyscale erosion. In simple EPS substitution, the original input multi-band grey levels 
of image 100 are simply replaced by the band EPS-filtered grey levels of smoothed image 142 wherever a 
pixel is marked as a defect in the classmap 175 (as shown in Figure 2). This method has the advantage of 
computational simplicity in that EPS values have already been computed during defect detection; however it 
requires that the EPS image be stored for use during cleaning. In addition, EPS substitution is relatively 
forgiving of "false alarms", i.e. pixels classified as dirt artifacts but which are really part of the scene. 
Conversely, EPS-substitution cleaning may not completely remove all of an extended defect, to the degree 
that the defect itself influenced the local EPS values. 

In the greyscale erosion cleaning option, an iterative process occurs over detected defect regions, 
whereby GLs of defect pixels are replaced by a weighted average of original GLs of "normal" (i.e. non- 
defect or corrected defect) pixels within a small distance of the defect pixel. As the iterations progress, the 
effect of the "normal" GLs propagates toward the center of the original extended defect. The iterations end 
when all marked defect pixels have been replaced. Two precautions are in order. First, the "normal "-labeled 
pixels that are in immediate contact with the defect pixels at the beginning of the iterations are also 
relabeled to class "defect" prior to the onset of cleaning (via a morphological dilation operation which uses 
a 3x3-pixel probe function). This is because such pixels, while failing the detection thresholds, often have 
some contribution from the defect in their GLs which will be visible when propagated over multiple pixels 
during erosion. Cleaning their values as well as those of "true" defect pixels can significantly improve the 
quality of the cleaning. The second precaution is that the spatial window over which the erosion weights are 
taken should be kept reasonably small, e.g. 5x5 samples. The reason for this is to minimize the influence 
of distant pixels which may not be part of the same scene structure as the region obscured by the central 
defect pixel . 

Current practice for the erosion cleaner is to use a 5x5-pixel sliding window with a minimum of six 
clean-class contained pixels in order for a central defect-class pixel to be filtered; if the central pixel is 
defect-class, but fewer than six clean samples are contained, then the defect classification is maintained to 
the next iteration. If the central pixel is of defect class and six or more clean-class samples are contained in 
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the window, then the GLs of the central pixel are replaced by a weighted average of the band GLs of the 
contained clean pixels, where the weights are linearly proportional to the inverse of the Euclidean spatial 
distance of the clean pixel from the central pixel. The cleaning is iterative (multiple passes through the class 
map) but is not recursive, i.e. the cleaning effect at one window location does not affect the cleaning at 

5 subsequent locations during the same iteration, in order to avoid directional biasing. 

Greyscale erosion has the advantage over EPS-substitution in being able to completely eliminate all 
trace of the artifact from the image; its disadvantages are increased complexity, a variable processing rate 
(the number of iterations is not known a priori), and the possibility for severe scene damage in the case of 
"false alarms" (i.e., real scene elements may be removed from the image). Despite these concerns, the 

w excellent cleaning performance of erosion is currently the preferred method in both Photo CD (Figure 7) 
and Cineon (Figure 8) applications. 

Alternative Implementations: 

75 The ADB algorithm has several "flavors". These include (a) operation at a single spatial resolution, (b) 
operation within a resolution pyramid, (c) single-defect-state detection, (d) multiple-defect-state detection, 
and cleaning via (e) EPS substitution or (f) grey scale erosion. In the single-state versions, all pixels of the 
image are classified as either "non-defect" or "defect", : i.e. only one defect state is treated. In multi-state 
versions, the "defect" classification is further divided into subclassifications depending on the confidence 

20 that a pixel is indeed a defect; thus a three-state classification would consist of "non-defect", "high- 
confidence defect", and "transition-confidence defect". In single-state versions the cleaner used is either 
EPS-substitution or grey-scale erosion, but not both in one implementation. The multi-state version may 
employ a combination of both EPS-substitution for cleaning the lower-confidence defect regions and grey- 
scale erosion for cleaning the higher-confidence defect regions. 

25 The flavor combining options (a), (c), and (f) is a simple variant to implement, though it is not 
necessarily the fastest or best-performing version. In this case no spatial resolution pyramid is formed, the 
residuals image is not segmented, and the cleaner is chosen to be gray scale erosion for maximum 
correction. Figure 7 illustrates such a process flow and is typical of a Photo CD application. In this case the 
preprocessing 110 consists of rescaling the pixel values of the original image for maximum viewing contrast 

30 by linearly remapping the lowest .5-% band gray level to a value of zero and the highest 99.5% level to 
255, as previously described in the description for Figure 2. The brightness test 130 for dirt on scanned 
negatives consists of minimum threshold pixel values which must be equaled or exceeded in each band of 
the feature image 120 in order for the corresponding pixel to pass the test as a possible defect. The 
residuals classification 150 is a single state decision per pixel as either a non-defect or defect class, based 

35 upon a requirement that the smallest band residual equal or exceed the zero-background residuals level 
and that the largest band residual equal or exceed some user-determined threshold level which is above the 
zero background level. The classmap conditioning 170 in this case consists of a morphological dilation of 
the defect pixels in the classmap using a 3 x 3 or 5 x 5 kernel. The cleaning 180 consists of greyscale 
erosion as previously described under the section called "cleaning. " 

40 Figure 8 illustrates an alternate processing flow which retains the characteristics, (a), (c) and (f), is 
appropriate for removal of bright defects in Cineon applications. In this case it is desirable to minimize the 
processing time required by the large image format the processing is therefor directed at only detecting the 
largest defects and at minimizing the number of bands in the feature image. The preprocessing 110 in this 
instance consists of creating a reduced resolution image rescaling it to a zero to 255 greylevel range such 

45 that the color balance is corrected, and then extracting a single band feature image in which the output 
pixel value is the larger of the Red, Green and Blue input values. The residuals test 150 applies a single 
threshold which each residuals pixel must equal or exceed to be marked as defect. Classmap conditioning 
170 consists of morphological dilation as described above. The cleaning module 180 consists initially of 
performing greyscale erosion of the reduced resolution unrescaled input image as previously described. 

so Next the cleaned reduced resolution multi-band image is bilinearly interpolated back to the original size and 
the conditioned classmap is replicated to the original image size. Finally those pixels in the original image 
which are marked "defect" in the replicated classmap are replaced by the pixel values in the interpolated 
cleaned reduced resolution image. 

It is possible to use ADB within a multi-resolution pyramid structure to provide some speed gains from 

55 use of multiple smaller EPS kernels. This implementation would also allow more effective cleaning from the 
use of EPS residuals thresholds which are tuned to the specific pyramid level. 

Figure 9 illustrates a flow for ADB processing within a resolution pyramid. In most respects, ADB 
cleaning within a spatial level is similar to that in a single-level ADB version, i.e. detection, classification, 
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and cleaning occur, and choices of one- or two-state and of EPS-substitution and/or greyscale erosion are 
the same. One key difference lies in the way the spatial pyramid residuals images are processed. 

In a spatial pyramid, an image is represented generally as a low-resolution image plus a number of 
residuals images of increasing resolution. To reconstruct the image to a particular resolution, one interpo- 

5 lates the low-resolution image by a prescribed method to the next resolution level and then adds the 
residuals values of that level to each pixel GL. This process continues until the desired resolution level has 
been reconstructed. The Photo CD file format is an example of a spatial pyramid structure. Note that these 
residual images are simply the pixel band difference between the original image at that resolution and the 
interpolated low-resolution version; they should not be confused with the EPS residual images discussed 

10 above. 

In a pyramid structure, all higher-resolution reconstructed pixels depend on the lower-resolution GLs at 
that spatial location. Thus when some subset of defect pixels is cleaned at a given resolution level, the 
residuals values at all higher resolutions at these locations must also be adjusted; failure to correct the 
residuals' values could result in "ghost" images of the original defects upon pyramid reconstruction. 

75 Although various residuals processing schemes are possible, the method currently used for ADB is the 
following relatively simple one: the base (lowest level) version of the original image is cleaned as in Figure 
7. The classmap 175 prior to base cleaning is then replicated in module 200 to the size of the next highest 
resolution level. The defect-class pixels in the replicated classmap are then dilated in module 210 by a 
spatial extent equal to the region of support of the pyramid smoothing function for that level, i.e., the half- 

20 width of the kernal used to create the next-lower-resolution level. The resulting dilated classmap 220 is then 
compared with the pyramid-residuals image 230; the pyramid residuals image is assigned a zero pixel 
value in every band wherever the corresponding pixel in the class map 220 is labeled as defect. The 
resulting modified pyramid-residuals image 250 is then added per-pixel in module 270 to a version of the 
base level cleaned image 190 which has been spatially interpolated to the present pyramid resolution level 

25 in module 260. 

This results in a reassembled (base + 1)-level image 100' in which the results of the base-level ADB 
learning have been incorporated. This image is then ADB-processed starting with the preprocessing step 
110', and the process of Figure 9 thus repeats until all levels have been reassembled and cleaned. Note 
that ADB cleaning may be omitted from higher-resolution levels at the expense of small defects not being 

30 corrected, but that all higher-level, pyramid-residuals must be zeroed as described in order to ensure that 
no defect "ghosts" will appear in the final image. 

While there has been shown what are considered to be the preferred embodiments of the invention, it 
will be manifest that many changes and modifications may be made therein without departing from the 
essential spirit of the invention. It is intended, therefore, in the annexed claims, to cover all such changes 

35 and modifications as may fall within the true scope of the invention. 
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Residual image 
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Spatial residuals zeroing block 
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Modified higher-resolution spatial resolution-residuals image 
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Cleaned-image interpolation block 
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Pixel value summing node 
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Rebuilt modified higher-resolution image 


110' 


Higher-resolution preprocessing 



The Invention may be summarized as follows: 

1 . A method for the detection and removal of defects in digital images comprising the steps of: 

a) preprocessing a digital image represented by pixel values to form a feature image in which defect 
pixels have large local contrast from their neighboring non-defect pixels; 

b) testing the value of each feature image pixel value to determine if the value is within a range of 
pixel values expected of defects; 

c) edge-preserving spatial filtering of each of the feature image pixel values that are within the range 
of expected defects; 

d) forming a residual value for each pixel as a function of the difference between the corresponding 
feature image pixel value and the filtered value of step c); 

e) testing each residual value of step d) to determine if the residual value is within a range of 
residuals values expected of defects; 

f) forming a map in which each pixel that is within the range of expected defects according to steps 
b) and e) is marked as a defect; and 

g) correcting the digital image using the map created in step f). 

2. A method for the detection and removal of defects in digital images comprising the steps of: 

a) preprocessing a digital image represented by pixel values to form a feature image in which defect 
pixels have large local contrast from their neighboring non-defect pixels; 

b) testing the value of each feature image pixel value to determine if the value is within a range of 
pixel values expected of defects; 

c) edge-preserving spatial filtering of each of the feature image pixel values that are within the range 
of expected defects; 

d) forming a residuals image which consists of the difference between the feature image pixel values 
and the filtered values of step c); 

e) testing each pixel value of the residuals image of step d) to determine if the residual value is within 
a range of residuals values expected of defects; 

f) forming a map in which each pixel that is within the range of expected defects according to steps 
b) and e) is marked as a defect; and 

g) correcting the digital image using the map created in step f). 

3. A method for the detection and removal of defects in digital images comprising the steps of: 

a) preprocessing a digital image represented by pixel values to form a feature image in which defect 
pixels have large local contrast from their neighboring non-defect pixels; 

b) testing the value of each feature image pixel value to determine if the value is within a range of 
pixel values expected of defects; 

c) edge-preserving spatial filtering of each of the feature image pixel values; 

d) forming a residual value for each pixel as a function of the difference between the corresponding 
feature image pixel value and the filtered value of step c); 

e) testing each residual value of step d) to determine if the residual value is within a range of 
residuals values expected of defects; 
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f) forming a map in which each pixel that is within the range of expected defects according to steps 
b) and e) is marked as a defect; and 

g) correcting the digital image using the map created in step f). 

4. The method wherein the map formed in step f) is conditioned by morphological dilation of the pixels 
marked as defect. 

5. The method wherein the map formed in step f) is conditioned by morphological erosion of the pixels 
not marked as defect. 

6. The method wherein the map of step f) is formed with each pixel marked according to the magnitude 
of its values according to steps b) and e). 

7. The method wherein the map formed in step f) is formed with each pixel marked according to the 
magnitude of its values according to steps b) and e)and is conditioned by morphological dilation of the 
pixels marked as defect. 

8. The method wherein the map formed in step f) is formed with each pixel marked according to the 
magnitude of its values according to steps b) and e) and is conditioned by morphological erosion of the 
pixels not marked as defect. 

9. The method wherein the filtering of step c) is performed with a k-nearest neighbor filter. 

10. A method for the detection and removal of defects in digital images represented by multiple spatial 
resolutions ordered from lower to higher, comprising the steps of: 

a) preprocessing a lower spatial resolution representation of the digital image by forming a feature 
image in which defect pixels have large local contrast from their neighboring non-defect pixels; 

b) testing the value of each feature image pixel value to determine if the value is within a range of 
pixel values expected of defects; 

• c) edge-preserving spatial filtering of each of the feature image pixel values that are within the range 
of expected defects; 

d) forming a residual value for each pixel as a function of the difference between the corresponding 
feature image pixel value and the filtered value of step c); 

e) testing each residual value of step d) to determine if the residual value is within a range of 
residuals values expected of defects; 

f) forming a map in which each pixel that is within the range of expected defects according to steps 
b) and e) is marked as a defect; 

g) correcting the lower resolution representation of the digital image using the map created in step f); 
and ' " 

h) correcting the higher resolution representations of the digital image using the map created in step 
0- 

1 1 . The method wherein one or more of the higher resolution representations previously corrected using 
the lower resolution defect map of step f) are again processed using steps a) through h) to detect and 
correct for defects too small to be detected at lower spatial resolution representations. 

12. The method wherein the filtering of step c) is performed with a k-nearest neighbor filter. 

13. The method wherein the filtering of step c) is performed with a k-nearest neighbor filter. 

14. An apparatus for the detection and removal of defects in digital images comprising: 

means for preprocessing a digital image represented by pixel values to form a feature image in 
which defect pixels have large local contrast from their neighboring non-defect pixels; 

first means for testing the value of each feature image pixel value to determine if the value is within 
a range of pixel values expected of defects; 

filter means for edge-preserving spatial filtering of each of the feature image pixel values; 

means for forming a residual value for each pixel as a function of the difference between the 
corresponding feature image pixel value and the filtered value from said filter means; 

second means for testing each residual value from said forming means to determine if the residual 
value is within a range of residuals values expected of defects; 

means for forming a map in which each pixel that is within the range of expected defects according 
to said first and said second testing means is marked as a defect; and 

means for correcting the digital image as a function of the formed map. 

15. An apparatus for the detection and removal of defects in digital images represented by multiple 
spatial resolutions ordered from lower to higher, comprising: 

means for preprocessing a tower spatial resolution representation of the digital image by forming a 
feature image in which defect pixels have large local contrast from their neighboring non-defect pixels; 

first testing means for testing the value of each feature image pixel value to determine if the value is 
within a range of pixel values expected of defects; 
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filter means for edge-preserving spatial filtering of each of the feature image pixel values that are 
within the range of expected defects; 

residual means for forming a residual value for each pixel as a function of the difference between the 
corresponding feature image pixel value and the filtered value of said filter means; 

second testing means for testing each residua! value from said residual means to determine if the 
residual value is within a range of residuals values expected of defects; 

mapping means for forming a map in which each pixel that is within the range of expected defects 
from said first and said second testing means is marked as a defect; 

means for correcting the lower resolution representation of the digital image using the map created 
by said mapping means; and 

means for correcting the higher resolution representations of the digital image using the map created 
by said mapping means. 

16. A method for the detection and removal of local defects in digital images comprising the steps of: 

a) generating a feature image from an original digital image; 

b) creating an EPS-residuals image using the feature image; 

c) testing and mapping image pixels based upon the values of the EPS residuals image and the pixel 
values of the feature image; 

d) modifying the mapped image as a function of mapped values of neighboring pixels; 

e) cleaning of defect pixels by an EPS-substitution process or by a greyscale erosion process; and 

f) cleaning of defect pixels in different spatial resolution representations of the original image if 
different spatial resolution representations are present. 
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Page 1 

Busteps 042193: (4-21-93) Applies a selective-averaging 

(E-NN smoothing) filter to a 3-band image; also • ■ , 

calculates the difference (residual image) between the original 
& the smoothed images - 

Search window size must be symmetric, odd of dimensions 

Filte^currently treats data as in range 0-4 999 *only*. 
Input image array is stored as Int*2. 

Intended for operation on a Sun workstation (e.g. Sparcstation 10) 
running Sun OS. 

Calls: 



75 



35 



40 



45 



50 



10 



Read Hdr isl 

"~ Timaae-f ormat -specif ic calls) 
Ihister * 

(image- format -specific caxls) 

Ihismnmx 
Winbust 042193 

— (image-format -specif ic calls) 

Hist Init 

Hist^Jpdate 

Iselavh 
Add Image 

~ (image- format -specific calls) 

program Busteps_0421 93 
integer*2 npix,nlin 

integer*2 ibuf f (4096 , 3) , buf f n (4096, 51 , 3 ) 
integer nav,npp, diagf Ig 

Integer npixm,nlm, np, nl,band, nbands, nbx 
Integer pixf onni, pixf ormo, pixf onnro, dumml / dumm2 
integer gmin (3) , gmax (3) 
integer hist (5000, 3) 
integer*4 britelo, britehi 

real*4 f actl , f act2, bias 

real*4 deltime, Dtirae, timearray (2) 

character*50 inf il, outf il, resf il 

data npm,nlm/51,51/ !Max EPS window size! 

Tdata pTxformro/3/ ! residuals image written as unsigned byte! 

data factl,fact2/l. 0,-1.0/ !Used in Add_Image ! 

data npixm, nbx/4096, 3/ 

data diagflg/0/ !For winbust 3! 

common /par /dununl , dumm2 , nav 



format ( 



Interactive parameter input 

write (6,5) 

"Busteps 042193":'/, ^ w . , . 

Tnput image may be 1,2, or 3 bands; /, 
Input values should be in range 0-4999; /, 
Max. image width: 4096 pixels per line; /, 
Wax. window height: 51 pixels'/, 
Left- & right -edge virtual mirroring /) 

write(6,10) ex 

f ormat ( 'Enter in filename I ch50] : ,5) 

read(5,' <a50) ' > infil !Must be ISL IPP imagefile! 

write (6, 30) 



55 
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JS 



35 



35 
40 



2 

VPS fllenamelchSO ]:'.$> 
formar( ' Enter out EPS tlx JSL ipp , 
30 read(5,' (a50)')outfil ,., S) 

wr!?e<6,35) resid uals filename ch50] . .5) 

sS2«* 4 s*.r .-.tf&sBs; sis--"' 

Call Exit(l) 
Endif ni=n? 

write (6,i?r nP * nl o of averag in 9 Poincs in window', 
fSrmatC Enter no. of aver g , m . n l)/2 ]:'.S> 

vrite(6,Sl> _ nl , ty pically about (np 

format \ 1 

stop 
Endif 



50 

formatC 
51 read<5,*> nav 



20 

60 
70 

25 



- ite '<h 6 E Lter low- . hi9h-bri9^n.« threshold, [intl 
Sfs.'MTrltelo.britehi 

r^rlnter out residuals bias constant:'.?, . ^ 

^'•r&isx.in.il.npi-n^-^^-^ 0 ^ 1 

^^~&>S£* — «• 1 

stop |5 

Endif ^ £ . 

. i • ~ a rravs & constants 
Initialize arrays * ^ 

30 npixm=npix 

„bx-nbands 9/g/92! 

deltine-Dtime(timearray) . 



pixfomo^ixfortni !Adci ed 8/T/ 321 

pixformro=pixformi 

, n Q f input image 
Calc histogram min,max oi i * 



""write (6,100) infil f 5Q) 

10 0 ft«tC inP-t *«• iS •»* ' -# ly fUt er the image. 

i-i Ac recruired to f r- 

Note: Tbis min/max call is r qu 

" lM ssSl1i- ko ' .«* ' a ""'' , 

A Do ib-l.nbands 
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20 



write (6, 151) ib, grain (ib) , gmax (ib) 
151 format (5x,i4, 2i9) 

Enddo 

..test if image mins & maxs are within LDT range of 0-4999: 
Do ib=l,nbands 

If (gmax(ib) .gt. 4999 J then 

write (6,*)' Warning — gmax .gt. 4999 STOP' 

Call ExitU) 
Elseif (gmin(ib) .It. 0) then 

write (6,*)' Warning — gmin -It. 0 STOP' 

Call Exit (1) 

Endif 
Enddo 

... Note: If gmax gt 4999, then either rescale input data, or redimension 
his arrays in selav3 & selavb, & reset threshold on above tester. 



Create EPS-filtered image 



write (6, *) ' [Beginning filtering] ' 
If (np . ge. 3 .and. np .le. npm) then 

Call Winbust 042193 <inf il , outf il, ibuff, buffn, 

1 npixm, nEx, nlm, np, nl, groin, gmax, pixf ormo, britelo, briteni, 

2 diagflg) 

Else 

Write (6,*) ' Np out of range stop' 

Endif 



Create EPS residuals image 

25 Call Add Iniage(infil,outfil f resfil,factl,fact2,bias,pixformro) 

Treslduals file is unsigned byte w/ mean^bias" ! 

write < 6, 200) outf il 

200 formate Ending; EPS file is ',a50) 
write<6,201) resfil 

201 formatC EPS residuals file is ',a50) 

30 

deltime=Dtime (timearray) 
write (6, 800) timearray (1) , timearray (2) 
80C format </' User time (sec):', 

x f8.2,'; system time (sec) : ' , f 8 .21) 

end 

35 



Ibismnmx: (5-28-91) Calcs. min i max values in a multiband histogram. 

Subroutine Ihismnmx ( hi st,nbands, gmin , gmax) 
40 Integer hist (5000, 3) , nbands, gmin (3) , gmax (3) 

Do ib«l, nbands 

Do k=l,5000 

If (hist()c,ib) .gt. 0) then 
gmin(ib)=k -1 
go to 50 

45 Endif 
Enddo 

50 Do x«5000,l,-l 

If (hist(k,ib) .gt. 0) then 
gmax (ib) «k -1 
go to 60 



Endif 

Enddo 

60 continue 
55 Enddo 
return 
end 
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job: read hdr isl.vf 

Date: Wed £pr 71 14:12:06 1993 
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25 



Read Hdr ISL: (9-25-90) Reads key info from header of an ISL-frat 
~ image header. Closes image f ile^upon^conclusxon^^^^^ 

Subroutine Read_Hdr_ISL <f ilename, npix, nlin,nbd,pxf orm) 

Integer imgnum, ios , npx , nln, nbd, pxf onn, rwmode 
IntegerM opnimg, getdef 
Integer*2 npix,nlin 
■Character*50 filename 

i — — — — 

include ' /local/include/iopacxagef ' 

rwmode =0 
imgnum^O 

ios=Opnimg (imgnum, filename, rwmode, . FAXSE . ) 
If (ios -NE. SYSNRM) then 

write(6 f *) ' Error during OPNIMG— Return' 

return 

Endif 

ios=Getdef (imgnum, npx, nln, ni>d,pxform) 
If (ios .NE. SYSNRM) then 

write (€,*)' Error durang Getdef — Return 
. return 

Endif 

Call Cls imq (imgnum) 
npix=npx 
nlin=nln 

return 
end 
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- Ihister: (5-28-91) Calcs. multiband histogram of an image using a line 
; buffer only. Intended for use w/ select ave2 .for . 

Subroutine Ihister (infil, npixm, nbx, ibuf f ,hist) 
Character*50 infil 

5 integer arrtyp?ostatus, gstatus, rstatus, imgin, npx, nln, nbands, pixf orm 

integer opnimg, get def, rdline, line 
integer hist (5000,3) 
integer*2 ibuf f (npixm, nbx) 

: include ' /local/include/iopackage.inc' 

w arrtyp-IDINT2 

-....Initialize the histogram: 
Do ib«l,nbx 

Do *=1,5000 

hist (x, ib)=0 

Enddo 

75 Enddo 

. . en the image: 
imgin*0 

ostatus=opnimg (imgin, infil, 0, .false . ) 
If (ostatus .ne. SYSNRM) then 

write (6, 15) infil , ^ , . 

20 15 formate Ihister: input image ',a50,' not open—abort ) 

stop 

Endif 

gstatus=*getdef (imgin, npx , nln, nbands, pixf orm) 
If (gstatus -ne. SYSNRM) then 

write(6,20) , ^ . 

20 formate Ihister: Gstatus error on input — abort ) 

stop 

Endif 

:...Read the image lines into the histogram: 
Do j=l,nln 

rstatus=rdline (imgin, line, -1, ibuf f, npixm, arrtyp) ■ 
If (rstatus -ne. SYSNRM) then 

write (6, 50) j \ 
50 formate Ihister: error reading line ,i4, —abort j 

stop 

Endif 

Do ib=l, nbands 

Do i=l,npx 

k=ibuff (i.ib) +1 
hist (x,ib)=hist <k,ib> +1 

Enddo 

Enddo 

Enddo 

Call clsimg (imgin) 
return 
end 
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Winbust 042193: _ . 

Uses a barrel-shifting line buffer in multiple bands. 
Implements a "fast-histogram 0 method. 



10 



25 



30 



Key Pa J^|£ erS £ h * 5 Q^ najne of input image file; not overwritten. 

outfil: Ch*50, name of output filtered image file- e ^ MH h« 

ibuff: lnt*2| all-band line buffer array; must be dimensioned by 
calling program. 
Int*2, strip buffer array; must be dimensioned by calling 

program. , , 

Int*4, Max. allowed width of image, m pixels per line. 
Input. , . . ^ 

Int*4. Actual width of image, m pixels per line. 
Read from image header. 

Int*4, Max. allowed no. of bands in input image. Ir.put.- 
Int*4, Max. allowed height of filter window (in lines) 
& thus max. allowed height of strip buffer. 



buffn: 

npixm: 

npix: 

nbx: 
nlm: 



np,nl: Int*4, 



Input . . . , „ 

Actual width, height of filter window in pixels, 



lines. Input. 

gmin,gmax:lnt*4, 3-eleraent vectors. Assumed mm & max band GLs 

y of input image. Assumed previously computed. Input. 

Note this limits current max. no. of bands to 3. 
pixformorlntM, Datatype of output filtered IPP image file. 

E.g. 3=uns byte, 4=signed Int*2, 6=signed R 4. 

Input. 

gmn,gmx: Int*4. Returned via Common. Set to current „ a1 „„ 

-band value of gmin & gmax, respectively. These values MUST 
be set prior to calling Iselavh; Function Iselavh uses 
these values passed via Common.. _ . _ _ . 

nav: Int*4, no. of samples to average m each filter window. 
Passed thru Common by caller. 

General comments: 

(1) 99% of the complexity of this subroutine is due to . 

the implementation of the strip buffer. The advantage of the 
strip buffer is that it allows "convolution "-type window 
operations while keeping at any one time m memory 
npix*nl*nbands worth of image, instead of npix*nlin*nbands for 
input & output images. rt ro« r.«c 

(2) Note that the EPS filter is *N0T* recursive, i.e. it only uses 
OLD? original GL values in its input weights; . thus . output from 
previous window filter do not affect neighboring window outputs 
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Subroutine Winbust 042193 (infil, outfil, ibuff , buff n, npixm, nbx, 
1 nlm, np, nl , grain, gmax, pixf ormo, britelo , brit ehi , diagf lg ) 

Character^SO infil, outfil 

integer imgin , imgout , ost atus , gst atus 

integer nbx, nbands, pixf orm, pixf ormo . vH 

integer pstatus , ost atus2 , rstatus , wstatus , npx , nln , npixm , npix , npixo 

integer nlin f nlinb 

integer arrtyp, lininc, jtype 

integer gmn, gmx, Iselavh, nav, gmin (3) , gmax (J) 

integer nlm, np, nl , nl2 , np2 , nbinx, gl_cent 

integer opnimg, getdef , putdef , rdline, wrlme 

integer diagf lg 

integer*4 britelo, britehi, count 

integer °2 ibuff (npixm, nbx) , buff n (npixm, nlm, nbx) 
integer*2 his (5000) 
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integers els (4096) !Note implicit max. of npix=4096 ! 

data nbinx/5000/ 
5 common /par/cmn,gmx / nav 

include ' /local/include/iopackage . inc' 
: remove above comment -out before compilation! 

n^nT/2 DINT2 !Half-height of filter window 
10 !S.g, fo r ni=5, nl2=2 

np2=np/2 ! Half -width 

count-0 



75 



"..O^en'Input & output images for sequential I/O: 
imgin-0 

ostatus=SYSNRM ... — * 

ostatus=opnimg (imgin, mril, u# • J - aise - i 
If (ostatus .ne. SYSNKM) then 

S«t^ 7 Si»bS"i042193:ln image '..50.' not open-abort' ) 
20 ' u call Exit (1) 

Endif 

MtttSs-getdSf (imgin, npx, nln, nbands , pixf orm) 
If (gstatus .ne. SYSNRM) then 

write (6, "75) ^^-r* \ 

formate Gstatus error on input aoort ) 

Call Exit(l) 
Endif 

noix-npx - ! Actual no. pixels/line in image 

If ,„pixt r . R - 6 r>pix.l i tt.; i: _ o<2l53! Bpliilj , 

' Call Exit (1) ' 
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Endif 

nlinb-nlin +2*nl2 

! nbands ^Actual no. of bands m image 
^tt/'°-Winbus t 042193" : Npix Hlin Type Nbands Typeo Np NaV 

fonaac c2£SSi!f SStjS^. IB, 

imgout=0 

pstatus=SYSTTRM . « 

pst atus«put def < imgout , npx, nln , nbands , pixf ormo ) 

?f (PS " tf r!t^>f' S S^ e exror for output-abort' 
Call Exit(l) 

Endif 

OStatUs2-SYSNRM 

ostatus2=opnimg (imgout , out f il, 1 , . false. ) 

i "i" for write-only: 

Tf /osratus2 .ne! SYSNRM) then 

If < osra ™ e - (6 ,* r output image not open-abort' 
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Endlf 



Call Exit (1) 



12 



lininc=nlin/8 IVariable used in writing current line-number 
jtype-lininc ! Ditto (just screen i/o related) 



rstatus=SYSNRM 
wstatus=SY5NRM 



W 

110 



= Load 1st strip into barrel buffer: 

Do i-l,nl2 +1 !Do per line in 1st buffer strip: 

line=j -1 ! (Damned -zero-starting" C array convention) 

rstatus=rdline (imgin, line, -1, ibuf f , npixm, arrtyp) 
If (rstatus .ne. SYSNTRM) then 



write (6, 110) j , . , , M . . - . -, . 

formate "Winbust_0421 93" : Head error at line i4, -abort ) 
Call Exit(l) 
Endif 

•Transfer GLs from line buffer to strip buffer: 

i (someone should be able speed up all this moving around of data) 
?5 !i =horiz pixel index for original image; , 

* i2=horiz pixel index for central (non-mirror) pixels an strip t>\ 
!i3=horiz pixel index for edge (mirror) pixels in strip buffer. 

j2=j +nl2 
j3=nl2 +2 -j 
Do ib=l,nbands 

20 Do i=l,np2 +1 !Left 

i2=i +np2 
i3=np2 +2 -i 

buf fn (12, j2, ib) =ibuf f (i, ib) 
buf fn (i3, j2, ib)=ibuf f (i, ib) 
buf f n £i2, j3, ib) =ibuf f (i, ib) 
buffn (13, j3, ib)=ibuff (i, ib) 

25 Enddo 

Do i=np2 +2,npix -np2 -1 ! Central 

i2-i +np2 

buffn (i2, j2, ib)=ibuf f (i, ib) 
buffn (i2, j3, ib)=ibuff (i, ib) 

Enddo 

. Do i=npix -np2,hpix ! Right 

30 i2=i +n P 2 

i3=2*npix +np2 -i 

buf fn (i2, j2, ib).=ibuff (i, ib) 

buffn (i3,j2,ib)=ibuff(i,ib) 

buffn (12, j3,ib)=ibuff<i,ib) 

buf fn {i3, j3, ib)*ibuf f (i, ib) 

Enddo 

35 Enddo 

•Write a status display to standard output! 

If ((j -nl2 -1) .ge. jtype) then eHio 
write (6, 600) j !Write current line number to stdio 

jtype=jtype +lininc 

Endif 

40 Enddo !End top a j° loop! 

jend=nl '"lend" is barrel index to most recent line 
1c=jend -nl2 l^jc" is barrel index to "middle* line 
Above sections only serves to fill the strip buffer for the 
! first time; it gets more complicated as the buffer contents move 

! "down" in the image. 

45 
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the buffer line index is rotated. 

Do J-nl2 + 2,nlin +n!2 .1 . ! Indexes the last ; line. + 1) of current 

5 u J j 0= j - n i2 -1 ! Image index of output lane. 

!Mark which pixels in the current line will be filtered: 
Do i=np2 +l,npixb -np2 

Do ib=l» nfcands . _ , _^„ A 

If (buffnU, jc,ib> -GT. britelo .and fiir er « 
w . buffn(i, jc,ib) .LT. britehx) then IDon t filter. 

A cls(i)«0 

go to 210 

Endi f 

Enddo 

els(i)=l Iraarlc for filtering! „_„^.-t 

countUount +1 -Increment to-be-filtered pixel counter! 
15 210 continue 
Enddo 

Do ib-l,nbands 

gmn»grain (ib) 
gmx=gmax ( ib ) 
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!Load window in central part of_ba rr el_buffer^ 

i Hot^a^Ix^ within 1/2-window width of vert image boundaries 
iBeSow ^"ame'as above -12". i.e. strip buffer horiz index. 

i Center-of -window index, 1st window! 
io=i -np2 ! Output pixel index! 

Call Hist init <nlm,np,nl,np2,buffn,npix2n,nbx, 
npixb,nbands, i, jend, ib,his, nbmx) 

gl _cent=buffn,i /j c,ib, ! GL ^wind oocenter 

" <ClS< ibuf! < ?ioHb?=Is e lavh(his,nbinx. g l_cent) -Filter 



Else 

ibuff <io, ib) «=gl — cent 

35 Endif 



•Update histogram for later window locations in a line: 

Do i=n?2 +2,n P ixb -np2 ! Center- of -window index! 
io .J fl n p2 ! Output pixel index! 

Call Hist update(nlm,np,nl,np2,buffn, nhinxl " 

npixH, nbx , npixb , nbands , i , ^end, ib, his , nbinx ) 

•If passes brightness thresholds test, then filter 
gl.cent-buffnX^ib, la^lado^cj^.r 

If <C ^uxx(io?ib)ilselavh^^^ 



Else 
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ibuff (io, ib)=gl cent 
Endif " 14 



Enddo 



Enddo I End band loop! 

line«=jo -1 

wstatus=wrline (iragout, line, -1, ibuf f , npixm, arrtyp) 
If (wstatus .ne. SYSNRM) then 
write (6, 410) jo 

formate *Winbust_042 1 93" : write error at line i4, ' -abort ' ) 
Call Exit (1) 

Endif 

♦Rotate barrel-buffer line indices: 

jend=jend +1 

If {jend .gt. nl) jend=l 

jc=jc +1 

If (jc .gt. nl) jc=l 

If (jo .le. nlin -nl2 -1) then 

j2=jend IFor bottom mirror! 

!Read new image line: 
line=j -1 

rstatus^rdline (imgin, line, -1 , ibuf f , npixm, arrtyp) 

If (rstatus .ne. SYSNRM) then 

write(6,430) j w 
formate "Winbust_0421 93 - : Error reading line i4, -abort 
Call Exit (1) 

Endif 

Do ib=l,nbands ! Overwrite oldest line in barrel-buffer" 
Do i=*l,np2 +1 'Left 
i2=i +np2 
i3=np2 +2 -i 

buffn (i2, jend, ib) -ibuf f (i, ib) 
buffn (i3, jend, ib)=ibuff (i, ib) 

Enddo 

Do i=np2 +2,npix -np2 -1 ! Central 

i2=i +np2 

buffn (i2, jend, ib) =ibuf f (i, ib) 

Enddo 

Do i=npix -np2,npix ! Right 

i2=i +np2 
i3=2*npix +np2 -i 
buffn (12 , jend, ib) =ibuf f ( i , ib) 
buffn (13, jend, ib) =ibuf f (i, ib) 

Enddo 
Enddo 

Else ! Virtual mirror of bottom of image: 

j2=jend +1 IBuffer index of "top" line in buffer! 

If <j2 .gt. nl) j2=l 

j2=j2 -1 ! (12-18-92) 

If <j2 .It. 1) j2-nl 1(12-18-92) 
Do ib=l,nbands 

Do i~l,npixb 

buffn <i, jend, ib) «buf f n (i, j2, ib) 

Enddo 
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€00 



800 



1A 

Enddo 
Endif 

.Write a status display to standard output! 
If (j -ge. jtype) then 

SEI&S™ Jin. -,i4.' is completed') 

jtype«jtype -rlininc 

Endif 

Enddo I End main • j- loop! 



-Call clsimg(imgin) 
Call clsimg (imgout) 
return 
end 
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17 

Hist: init: (8-12-91) Used in WinbustJ)42l93. Initializes * calcs a histogram 
" of a window. 

Subroutine Hist_init (nl*, np, J 1 / n P 2 f£V^^ n P ixm ' nbX ' 
! npix,nbands, i, jend, ib, his, nbinx) 

integer*2 buf fh (npixm, nlm,nbx) 
integer*2 his (nbinx) 

integer nlm, np, nl, n P 2 , npixm, npix, nbx f nbands, i , jend, ib f nbinx 
integer gmn,gmx, jO 

common /par/gmn, gmx 

...Initialize histogram to zero: 
Do Jc=gmn+1, gmx+1 
his(k)=0 

Enddo 

...Load histogram of window: 
jO=jend 
Do jj«l,nl- 

10=j0 +1 

If (10 .gt. nl) j0=l 
Do il~l,np * 

i0=i -np2 -1 +ii 

ia~buffn<iO, jO,ib) +1 

fcis(ia)~his<ia) +1 

Enddo 

Enddo 

return 
end 
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Job: hist update -vf 

Date: Wed Spr 21 14:30:52 1993 
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integers buf fn (npiwn, nlm, nbx) 
integers his (nbinx) 



»*T&ft•** , * , * , ■ 1,, * U " ,,1 " 

integer jO, idrop, iadd, npix 



j 0 =jend 

idrop^i -np2 -1 
iadd*=i +np2 

Do Ij^lfHl 

ia-buffnUdrop, jO,ib) +1 
his(ia)«his(ia) -1 
ia*buffn(iadd, jO.ib) +1 
his(ia)«fcis(ia) +1 

Enddo 

return 
end 
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Iselavh: (6-19-91) Mod'd from selavb to eliminate heavy central 21 
weighting on initial-point in interest of speed. 
Should be valid for all window sizes. 

Note that this function runs more slowly as the image dynamic range 
[graax -gmin] increases. 
Assumes histogram is precalculated. 

°ac" is the GL of the "center" pixel in the current window. 
References: 

(1) °Noise Cleaning by Iterated Local Averaging", by L. Davis 
6 A. Rosenfeld, IEEE Transactions Systems, Man, Cybernetics, 
Vol. SMC-8, Wo. 9, pp. 705-710, Sept. 197B. 

(2) " Digital image smoothing and the sigma filter", 

by J. S. Lee, CVGIP, V. 24, pp255-269 (1983). Sigma filter, however, 
uses a ft fixed* range, in terms of +/- a*sigma, where sigma=est. 
std. dev. of additive noise. 

This function computes a selective-averaging filter output for a 
current location of a rectangular window of 1 band of an image array, 
S returns this output as the standard function return. 

Key Parameters: 

his: Int*2 (nbinx) ; Histogram of GLs in current window; 

Assumes that GL value corresponding to a histo bin is 
offset by 1, e.g. his(l) corresponds to GL=0. 
20 nbinx: Int*4; Max. no. of bins in histogram. 

ac: Int*4; spatially-central GL in current window. 

gran,gmx:IntM; min & max GLs in entire image for present band. 

Passed from caller via Common, 
nav: Int*4; desired no. of pixels to average in the 
current window. 

25 Output "iselavh" will be a weighted sum of histogram bin occupancies; 

this weighted sum is num/den. 



30 



Integer Function Iselavh (his, nbinx, ac) 

real* 4 val, vall,val2 

integer*2 his (nbinx) 
integer* 2 cnt . 

integer groin, gmax, nav, a c, nbinx 
Integer ar, al,denr, denl 
INTEGER* 4 NUM,DEN 

35 common /par /gmin, gmax, nav 

common /par2/num, den,ar, al, cnt, vail, val2,val 



S average the nav-nearest grey levels — 
den=his (ac +1) 

if (den .ge. nav) then • !If have >den local vals =ac,then no-change! 
40 iselavh=ac 
return 

endif 

num^ac^his (ac +1) 
do k°l,4999 

ar°ac +k ! Current GL to "right" of center GL 

45 al°ac -k ! w " -left" of - 

if (ar .le. gmax) then 

num»num +ar*his (ar +1) 
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endif 
if (al 



den=den +his (ar +1) 

. ge. gmin) then 
num=num +al*his (al +1) 
den =den +his (al +1) 
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endif 
if (den 



.gt. nav) then 



cnt=den 
denr=0 
denl=0 
If (ar 
If (al 



-nav 



else if 



le. gmax) denr=his(ar +1) 
ge. gmin) denl=his(al +1) 
va!2«Float (denr +denl) 
vall=Float (ar*denr) +Float <al*denl) 
val=vall/val2 

iselavh^Nint ( (Float (num) -cnt*val) /Float (nav) ) 
return 

(den ,eq„ nav) then 

iselavh=Nint (float (num) /float (den) ) 
return 



enddo 



Endif 

write (6, 
write (6, 
write (6, 
write (6, 
write (6, 



Call Exit (1) 



)' Unexpected termination in iselavh!!' 

)' Local Iselavh values at Exit:' 
) ' Central GL value (ac) : ' f ac 

)' Gmin, Gmax :', gmin , gmax 

) ' Nav: ' , nav 



end 
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24 

Add_Image: (6-26-91) Adds 2 1-3 band images w/ additive 
constants. 

Assumes the internal image representation is Real*4 
Max. 0 pixels/line-4096. 

Subroutine Add Image (inf ill, infil2, outfil, factl, fact2, const, pixformo) 
Character*50 inf ill , inf il2, outfil 
Real*4 rbuf f 1 (4096, 3) , rbuf f2 {4096, 3) 
RealM factl ,fact2, const 

integer imginl, imgin2, im^out, nimx, ostatus, gstatus 
integer nbands , nbandsc , pixf orm, pixf ormc, pixformo, band 
integer pstatus, ostatus2 , rstatus, wstatus, npx # nln, npx2 , nln2 
integer opnimg, getdef , putdef , rdline, wrline 
integer arrtyp, lininc, jtype, imax 

data imax/4096/ 

include ' /local/include/iopackagef ' 
arrtyp«IDREA4 

-..Open input color Image: 
imginl=0 

- ostatus=opnimg (imginl, inf ill, 0, .false. ) 
If (ostatus -ne. SYSNRM) then 
write (6, 10) infill 
10 formate In color image ',a50, ' not open-abort') 

Call Exit(l) 

Endif 

gstatus=getdef (imginl, npx, nln, nbands, pixf orm) 
If (gstatus -ne. SYSNRM) then 
write(6,20) 

20 formate Gstatus error on input — abort') 

Call Exit(l) 

Endif 

If (npx .gt . imax) then 

write (6,**) ' Npx .ge. imax — abort' 

write (6,*)' Npx=',npx 

write ( 6 , * ) * lmax= * , imax 

Call ExitU) 
Endif 
npix=npx 
nlin=nln 
write (6, 45) 

45 format (/' "Add Image": Npix Nlin Type Nbands Typeo' ) 

write {T, 50) npix, nlin, pixf orm, nbands, pixformo 
50 format (13x,2i5,i6, i8,i8> 

...Open 2nd input class map: 
imgin2°0 

ostatus^opnimg <imgin2, inf il2, 0, .false. ) 
If (ostatus .ne. SYSNRM) then 

write (6, 60) infi!2 
60 format e In 2nd image ',a50,' not open-abort') 

Call Exit(l) 

Endif 

gstatus«getdef (imgin2, npx2, nln2, nbandsc, pixf ormc) 
If (gstatus .ne. SYSNRM) then 
write (6, 65) 

65 formate Gstatus error on 2nd input — abort') 

Call Exit(l) 
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"Iki (#> ' "EUS'i-SifdiS^SiSl "un^al-CaU Exit: (1, • 

Call Exit (1) 
Endif 

P status4uldef Umgout, npx, nln. nbands, pixf ormo) 

Call Exit(l) 

Xtitus2-opni»g Umgout ^outf il^l alse . ) 

If (ostatus2 .nel SYSNKM) then aK rtrr , 

ir < OSM ™5 te(6# * r output image not open-abort 

Call Exit (1) 

Endif 

write(6,*>' Facrl. Fcat2, Const*' , factl. fact2, const 

Read each line of input image: 
lininc=nlin/4 
jtype=lininc 
Do j=l,nlin 

line=j-l 

•Read multiband line from c0 } 0 ***2f arrtV p) 
rstatus-rdline (imginl, line, -1, rbuf f 1, imax, arrtyp) 

i for all bands! 

If (r status .ne. SYSNRM) then 

?o^tr 8 Error reading image line -abort') 
Call Exit(l) 
Endif 

!Read line from l-kand c lassma £: arrtvP) 
rstatus-rdline (imgin2, line, -1, rbuf f 2, imax, arrtyp) 

t "-r for all bands! 
If (rstatus .ne. SYSNRM) then 

format r 8 Irror reading image line -abort') 
Call Exit(l) 
Endif 



Do band-1 , nbands 
Enddo 

...Clip data to output range: 
If (pixf ormo .eq. IDBYTE) then 

D ° iel ' n &ufflU,band^ } 
rbuff 1 (i,band) -Aminl (rbuf f 1 (i,band> ,255.) 

Elsl n if° (pixf ormo .eq. IDINT1) then 

i — 1 rir-> 4 v 



45 



50 



Do i=l , nplx band) ffl (± , band) , -128 ) 

rbuf f 1 <i, band) -Aminl (rbuf f 1 (i , band) , 127 . ) 

Else^if 0 (Pixf ormo .eq. IDINT2) then 

Do iral ' n ^ ff i (irba nd)^int(rbuffia,^ 

rbuffl(i,band)»Amaxl(rbufflU, band) ,-32768.) 
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70 



rbuffl (i,band)=Aminl (rbuffl <i, band) , 32767 . ) 26 

Enddo 

Else if (pixformo .eq. IDREA4) then 
Else 

write (6, 100) s , . , \ 

forraatC Pixf ormo=' , 13, ' not implemented-abort ) 
Call Exit(l) 
Endif 

Enddo !End band loop! 



c ...Write output line of data: 

wstatus^Wrline (imgout, line, -1, rbuffl , imax, arrtyp) 
If (wstatus -ne. SYSNRM) then 

write<6,*)' Error in vstatus-abort' 
. Call Exit (1) 
Endif 

If (j -ge. jtype) then 

write<6,200) j t 
'5 2G0 format (I0x,* Line ',14,' is completed ) 

jtype= jtype +lininc 

Endif 

Enddo !End line loop! 

c Close images 

Call cisimg(imginl) 
20 Call clsimg(imgin2) 

Call* els img (imgout) 
return 
end 
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Appendix B-Dockfet 66,719 



© 1993 Eastman Kodak Company. 
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c Bust Class Creates ** two-state ** class map from 3-band residuals image, 

c "Optionally dilates classmap before output. ,m ™ t * - * • i 

c Calls: Irgb_load_isl, Buster_Class_Create; Dilate_Class, Ioutput_isl. 

c Image arrays: . ' 

c ' res : 3-band dirty residuals image (input) 

c els :i-band classmap (output) 

c residuals image & class map are assumed lnteger*2 values (i.e., not real 

c Note demonstration max image size of 1B70 x 1374. 

c Does not assume that residuals image is biased by a predetermined 

c constant; this info is implicit in the GL-res thresholds entered by the 



c Calls: 

c Irgb load isl . . • % 

c ~ (calls image-format-specific 10 routines) 

c Buster_Class_Create 

c Dilate Class 

c I output isl t , . 

c "(calls image-format-specific io routines)^ 



20 



Program Bust__Class 

Integer*2 res (3, 1870, 1374) 

Integer*2 els (1870, 1374) , cls2 (1870, 1374) 

integer*2 npix,nlin 

Integer npixm, nlinm, band, nbands , nbc , nbx , cnt , np , ni 
integer pixf ormi,pixf ormo 

integer cnt tot 
integer ixO, iy0,nx, ny 

integer dgl , glmin , glmax , gltmm , gltmax ( dilf lg , byte 

real°4 deltime,Dtime, timearray (2) 

character*50 inf il, outf il, resf il, clsf il 
character* 1 iname (50) ,iname2 (50) 
character*50 commSO 

equivalence (comm50, iname2 (1) ) 

data dgl/255/ 'Defines "dirt" class! 

data byte/ 3/ 

data nbx, npixm, nlinm/ 3, 1870, 1374/ 

data pixfonno/3/ ! Classmap of type "byte 

i — — — — — ~ ~ 

write (6,°)' (Max. image si2e: 1870 x 1374]' 
write (6,*)' ' 

write (6, 20) , r„v.cm.' *\ 

formate Enter in biased residuals image filename [ch50] . 
read (5,' (a50) ' ) resfil 



write (6, 30) j , ei 

30 formate Enter out classmap image name: ,5) 

read<5,' (a50) ' ) clsfil 

write (€,50) t , . , e , 

50 formate Enter min, max residuals GL thresholds: ,sj 

read(5,°) gltmin, gltmax 

write (6, 60) 
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60 formate Do you wish prefilter dilation? Eyes="l"] : ' , 

read (5,*) dilflg 

write (6,*) ' Class threshold used is ',dgi 
deltime=Dtime (timearray) 

■Read biased color-residuals image: 

Call Irgb_load isl (resfil , res, rtbx, nbands, npixm, nlinm, 
1 npix, nTin, pixformi , 0) 

write <6, 150) resfil 
150 formate Input residuals file is ',a50) 

i _ — 

! Compute class map: 

Call Buster_Class_CreaLe (res, cis, nbx, nbands, 
1 npixm, nlinm, npix, nlin, cltmin, gltmax, dgl) 

i 

■Optionally, dilate classmap by one pixel: 

If (dilflg .eq. 1) then 

write (6,*)' [Dilating class map]' 

Call Dilate_Class (els, cls2, 
1 npixm, nlinm, npix, nlin, dgl) 

Call Ioutput__isl (clsf il, cls2, npixm, nlinm, 
1 npix, nlin, byte, iname) 

write (6,200) clsfil 
200 formate. Out dilated class file is ',a50). 

Else ! write undilated classmap! 

Call Ioutput_isl (clsfil, els, npixm, nlinm, 
1 npix, nlin, byte, iname) 

write (6, 300) clsfil 
300 formate Out undilated class file is ',a50) 

Endif 

deltime=Dtime (timearray) 
write (6, 800) timearray ( 1 ), timearray (2) 
800 format (/' Bust Class:: user time (sec):', 
1 f8.2,'; system time (sec) : ' r f 8.2/) 

end 
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c Irgb load_isl: Loads up to a 3-band image into a lnt*2 array. 

Subroutine Irgb_load_isl <bname , b, nbx , nb, npixm, nlinm, npix , nlin , pixf orm, 
1 lunprt) 

Character*50 . bname, inf il ^ v 
Integer*2 b (nbx, npixm, nlinm) ,npix, nlin, ibuf f (4096,3) 
integer imgin, gstatus, ostatus 
integer nbands, pixf orm, band, nbx, nb, npixm, nlinm 
integer rstatus , npx , nln 
integer opnimg, getdef , rdline 
integer arrtyp, lininc , jtype , pixdim 

Data pixdim/4096/ 
include ' /local/include/iopackagef ' 

arrtyp=iDlWT2 

c Open input image: 

15 imgin^O 

inf il=bname 

ostatus«=opnirog (imgin, inf il, 0, . false . ) 
If (ostatus .ne. SYSNRM) then 
write (6, 10) infil 
10 formate In image ',a50 f ' not open-abort' ) 

stop 

20 Endif 

gstatus^getdef (imgin, npx, nln, nbands, pixf orm) 
If (gstatus .ne. SYSNRM) then 
write (6, 20) 

20 formate Gstatus error on input — abort ) 

stop 

25 Endif % k ^ 

If (nbands .gt . 3 .or. nbands .gt- nbx) then 

write (6, Input no. bands too large — abort' 
stop 

Endif . ^ 

If (npx .gt. npixm .or. npx .gt. pixdim) then 

write<6,*)' w IRGB_load_isl w : too many pixels — abort' 
30 stop 
Endif 

npix=npx 
nlin=nln 
nb=nbands 

write (6,45) npix, nlin, nbands, pixf orm „_,••» 
35 45 format(2x,' "Irgb_load_isl" : : Npix:',i5,'; Nlin:',i5,'; Nb: ,i2, 
1 '; Type:',i2) 

If (lunprt ,gt. 0) then 

write (lunprt, 45) npix, nlin, nbands, pixf orm 

Endif 

c Read each line of input image: 

40 lininc=nlin/4 
jtype°lininc 
Do j«l,nlin 

c ... Read input lines & load 3-d array: 

line=j-l 

rstatus=rdline (imgin, line, -1 , ibuf f , pixdim, arrtyp) 
! for all bands! 

45 If (rstatus .ne. SYSNRM) then 

write (7, 80) j 

80 formate Error reading line -abort') 

Stop 
Endif 
Do i«l,npix 

50 
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Do band=l,nbands 5 
b <band, i, j ) =ibuf f (i, band) 

Enddo 
Enddo 

If <j .ge. jtype) then 
5 write (6,200) j 

200 format dOx,' Line ',i4,' has been read') 

jtype= jtype +lininc 

Endif 

Enddo !End line loop! 

c Close image: 

'0 call clsimg(imgin) 

return 
end 
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Job: buster class create dugas.vf 
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of defects. 

Subroutine Buster ClassJTreate (res. els, nbx, nbands, npixm, nlinm, 
npix, nlinTgltmin, gltmax, dgl ) 

Integer *2 res (nbx, npixm, nlinm) , els (npixm, nlinm) 
Integer*2 npix,nlin 

integer nbx, nbands, npixm, nlinm, gitalin, gltmax, dgl, glmin, glmax 



c 
c 
c 
c 



10 
11 



20 



25 



formal r^ust ClassJTreate: Gltmin Gltmax Dgl') 
write (6, 11) aTtmin, gltmax, dgl 
format ( 1 9x, 2i7 , i5 / ) 

Do j«l,nlin 

Do i~l,npix t 

glmin°res (1, i, j) 
g lmax^g lmin 
Do ib»2, nbands 

glmin°MinO (glmin, res (ib, i, 3J ) 
glmax»MaxO (glmax, res (ib, i, 3) ) 

?? d ?glmin .ge. gltmin .and. glmax .ge. gltmax) then 

cls(i, j)=dgl 

Else , n 

cls(i, j)-0 

Endif 

Enddo 

Enddo 

return 

end 
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c D i- t e s CIass: a Dilates f f ^s^indo*. 

c ctar-class which is dilated; if current pixel is undilated then 

c old classvalue is copied as new classvalue. 

c Also -Si^^S^Si^ still covering i»age. 

^rout i^^a^C^ nlinm, npix, nlin, ctar ) 



c 



,o Integers clsi (npixm p nlinm) , clso (npixm, nlinm) , npix, nlin 

Integer npixm, nlinm, ctar, im, ip, jm, jp 

i 

write <6, 10) ctar „ . , icy, 

format </' Dilate_Class: target class is ,15/) 



10 

'5 Do j*=l,nlin 

1m°Max0 (jm, 1) 
jp~j+l 

jp«MinO (jp,nlxn) 
Do i=l,npix 

im=i-l 

20 im*MaxO (im,l) 

ip=i+l 

ip=MinO <ip,npix) 
Do jj«1m, jp 
Do ii»im,ip 

!If any neighbor^ctar, then 
25 i set center^ctar as well: 

If (clsi(ii,jj) -eq. ctar) then 
clso(i, j)«ctar 
» 100 



30 



35 



40 



go to 
Endif 

Enddo 

Enddo _ ^ 

!No change, so copy to output 
clsoU, j)»clsi<i, j) 
100 continue 



Enddo 



return 
end 
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ll 

r SUBROUTINE: IOtn.*-UT isl (10-9-90) 

C PROGRAMMER: ELIZABETH J. DONALDSON/R. Gray 

C DATE: 9/20/90 

r OFILE - INPUT - CHAR*50 - DESTINATION IMAGE 

r FOUT - INPUT - Integer*2 - DESTINATION ARRAY 

I NPIXM - INPUT - INTEGER* 4 r MAX PIXEL DIMENSIONS FOR FIN 

C NLINM - INPUT - INTEGER* 4 - MAX LINE PENSION FOR FIN 

C NPIX OUTPUT - INTEGER* 4 - NUMBER PIXELS IN INFILE 

C KLIM - OUTPUT - INTEGER* 4 - NUMBER LINES IN INFILE 

r PIXFORM - OUTPUT - INTEGER* 4 - INFILE DATATYPE 

C INAME50 - OUTPUT - CHARACTER - (UNUSED) = 



c= 



= Subroutine I output ^isl (of ile, f out , npixm, nlinm, npix, nlin, pixform, 
+ inameSO) 

include ' /local/include/iopackage . inc' 

character of ile*50, inameSO (50) 

15 lnteger*2 fout (npixm, nlinm) 

integer*2 npix, nlin 

integer putdef , opnimg, wrband t 

integerM imgout , arrtyp, gstatus, pixform, npixm, nlinm 

integer*4 ostatus, pstatus, wstatus , npx^nln 

write (6,5) npix, nlin, pixf onn , , ± 

20 5 formatC I output isl: image header data: np « , H/ n± 

' pixform = ' , i3T 

arrtyp « IDINT2 
imgout =0 

C DEFINE CHARACTERISTICS OF OUTPUT IMAGE 

25 npx=npix 

nln=nlin . 
©status = putdef (imgout, npx,n In, 1, pixform) 
if ( gstatus .ne. SYSNRM) go to 1000 

C OPEN OUTPUT IMAGE FOR RANDOM ACCESS 

ostatus =» opnimg (imgout, of ile, 1, .true. ) 
30 if ( ostatus -ne- SYSNRM) then 

write (6, 10) of ile 
in format (' Output image' , a50, ' not open. ) 

write (6,*) ' Exiting — ostatus= ' , ostatus 
go to 1000 

endif 

35 - if ( pixform. eq.IDBYTE) then 

do 100 j = l,nlin 

do 200 i - l,npix 

if (foutti, j) -lt.0) fout(i,j) « 0 
if (fout <i, j) -gt.255) fout(i,j) = 255 
200 continue 
100 continue 
40 elseif( pixform. eq.IDINTl) then 

do 110 j * l,nlin 

do 210 i - l,npix 

if (fout(i, j) .It. -128) fout <i, j) - ;128 
if (fout <i,j) -gt. 127) fout U, 3) - 127 
210 continue 
110 continue 
45 elseif ( pixforra.eq.IDINT2) then 

do 120 j « l,nlln 

do 220 i ° 1 , np i x . _ rt 
if ( fout (1,1) .It. -32768) fout (i, j) ";32768 
if ( fout <i,j) .gt. 32767) fout (i,j) ~32767 
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220 continue 
120 continue 

elseif ( pixforro.eq.IDREA4) then 
GO TO 900 

else 

write (6, 30)pixform 
30 formate IOUTPUT_ISL: Pixform ',i3, ' not implemented.') 

endif 

900 ustatus » wrbandiimgout, 0, £out,npixm, nlinm, arrtyp) 
if ( rstatus .ne. SYSNRM) go to 1000 

1000 continue 

call clsimg (imgout) 

return 

end 



12 



15 



20 



Appendix C — Docket 66,719 
1993 Eastman Kodak Company, Rochester, W.Y. 14650-2201 
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Buster Erode am: (11-11-92) Uses dynamic memory allocation on the Sun. 
Calls HedirtSxS byte2 nopad; no 2-pixel -dead jone" 
around image; aTl pixels in image are potentially cleaned. 

Intended for use on a Sun workstation running SunOS (e.g. a 
Sparcstation 10) . 



5x5 erosion window is currently hard-wired in. Minimum clean fraction 
of erosion window is nominally .24. 

Cleans RGB image by grayscale erosion 
using an external 2-state class map. 

Calls: 

Irgb load isl . „ „ s 

~~ (includes image-format -specific calls) 

Binput Isl 

~ (includes image-f ormat-specif ic calls) 

Dedirt5x5Jbyte2_nopad 

Irgb write isl 

~" (includes image-f ormat-specif ic calls) . 

Image arrays ; .... 

1ml zOriginal (dirty) 3-band image (input) 

els :l-band classmap (input) 

jm? : Cleaned 3-band image (output) 

Warning: Once erosion iteration begins, imi & im2 contents 
are ^partially* -cleaned image; original is not retained. 

Removes by interpolation all pixels whose class -map GL is .ge. 
a threshold value. Works on a multiband image (up to 3 bands) . 
Assumes iterative in/of recursive filtering. 

Uses adjustable variable-input window size, allows only filtering 

a rectangular piece of the image. 
Image £ class map are assumed Integer* 2 values (i.e^, not real)^ 



Program Buster_Erode_Am 

Pointer <piml # imi) , (pim2 ( im2) , (pels, els) , (pcls2, cls2) 

Integer*2 imi (1, 1, 1) f im2 (1, 1, 1) 
Byte cls(l,l), cls2(l,l) 

Integer*2 imi (3 f 1828, 1332) , im2 (3, 1828, 1332) 
Byte els (1828, 1332) , cls2 (1828,1332) 

integer*2 npix,nlin 

Integer npixm, nlinm, band, nbands, nbc, nbx # cnt 
integer pixf ormi, pixf ormo, pixf ormc 

integer iml_bytes, im2_bytes, cls_bytes, cls2_bytes, tot_bytes 
integer itmax, it, cnttot 
integer inflg 
integer dgl,byte 
Integer*4 iystart, iyend 

real 0 4 deltime, Dtime, timearray (2) , fmin 

character*50 inf il, outf il,clsfil 

data itmax/40/ !Max. no. of erosion iterations! 

data dgl/255/ !Defines "dirt* class! 

data dlagflg/O/ 
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data byte/3/ 

data nbx, npixm, nlinm/3, 1826, 1332/ 



write(6,*) # — — — — — — ™r 

write (6^*)' Buster_Erode_Am : Dynamically allocated memory' 

write (6, 10) , , 

10 format ( /' Enter in ^dirty" file name[ch50j: ,S) 

read (5,' <a50)')infil 

Call Read hdr isl (inf il, npix, nlin, nbands, pixf ormi) 
write (6, 15) npTx, nlin, pixformi, nbands 
15 formate npix=',i4,'; nlin«',i4,'; type=',i2,'; bands= ,12) 

npixm=npix 
nlinm=nlin 
nbx=nbands 

iml byt e s* 2 * nbx* npixm *nl iron 
im2~ljyt e s= iml_byt e s 
c 1 s~by t e s«np ixm * nl inm 

cls7 bytes=cls bytes m v , . . ^ 

tot Sytes«iml 5ytes +im2_bytes +clsjbytes +cls2 bytes 
write (6,*) '"Memory allocation requested (bytes): ',tot_bytes 

write<6,20) # ^ % 

20 formate Enter in class filename:',*) 

read (5,' (a50)') clsfil 

write (6, 30) r WCA , , e% 

30 formate ' Enter out "cleaned" filename [cnSOj : ,5) 

read<5,' (aSO) ' ) outf il 

write (6, 40) . ojC , „ , _ , c . 

40 formate Enter reqd cleaning fraction "fmin [Real] t ,>) 

read (5,*) fmin 

If (fmin .le. 0. .or. fmin .gt. 0.5) then 

write (6,*)' Fmin out of range — abort' 
Call Exit(l) 

Endif 

c. ............. 

write(6,*) ' Class threshold used is ',dgl 

!Read color image to be filtered: 

write (6, 150) inf il 
150 format </' Input (dirty) file is ',a50) 

deltime*Dtime (timearray) 

piml«malloc (iml bytes) 

If (piml .eq. 0J then w _ 

write (6,*)' Malloc for "iml° failed— abort 
Call Exit(l) 

Endif 

Call Irgb_load isl (inf il, iml, nbx, nbands ,npixxn,nl inm, 
1 npix, nTin, pixformi, 0) 

pixformo~pixformi !Out datatype set equal to in type! 

vrite (6,160) clsfil 
160 formate Input class file is *,a50/) 
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pcls*=raalloc (els bytes) 
If (pels .eq. OT then 

write <€,*)' Malloc for "els" failed— abort' 

Call Exit(l) 

Endif 

Call Binput isl (clsfil, els, npixm, nlinm, npix, nlin, pixformc, 
l,nEc) 

pim2**malloc (ixn2 bytes) 
If (pim2 .eq. OT then 

write (6, °)' Malloc for "intf" failed— abort' 

Call Exit (1) 

Endif 

pcls2==malloc (cls2_bytes) 
If (pcls2 .eq. 0) then 

write (6,*)' Malloc for a cls2" failed— abort ' 

Call Exit (1) 

Endif 

i — — 

inflg=l ! Start by assuming w cls a is input class map! 
! Clean by erosion iteratively until no more pixels are replaced: 
write (6,*)' [Beginning cleaning]' 
cnttot^O 

iystart^l IRequired initializations 

iyend=nlin 

Do it«l,itmax 
cnt=0 

If (inflg .eq. 1) then 

Call Dedirt5x5 byte2_nopad (irol, els, im2, 

cls2, npixm, nlinm, npix, nlin, nbx, nbands, 

dgl, cnt, fmin) 
Call Dedirt5x5 byte2_nopad2 (iml, els, im2, 

1 cls7, npixm, nlinm, npix, nlin, nbx, nbands, 

2 dgl, cnt , fmin, iy start, iyend) 
inflg«2 

Else if (inflg .eq. 2) then 

Call Dedirt5x5_byte2_nopad(im2,cls2, iml, 

els , npixm, nlinm, npix, nlin, nbx, nbands, 

dgl, cnt, fmin) 
Call Dedirt5x5_byte2_nopad2 (im2, cls2, iml, 

1 cls r npixm, nlinm, npix, nlin, nbx, nbands, 

2 dgl, cnt, fmin, iystart, iyend) 
inflg=l 

Else 

write (6,*)' Inflg out of range — abort' 
Call Exit(l) 

Endif 

cnttot=cnttot +cnt 
write (6,200) it, cnt 

formate After it=',i3,' pels cleaned = ',i6) 

If (cnt .eq. 0) then !No more dirt-class pels to be cleaned! 
go to 500 !Dreaded "go-to" statement! 

Endif 

Enddo a m , 

write(6,*)' Warning-Cleaning loop exhausted after it=' , itmax 
write(6,°)' Aborting due to loop exhaustion' 
Call Exit(l) 
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continue 

write <6, 550) it , . 

formate Cleaning ended after iteration-' , i3) 
write (6, 552) cnttot 

fonaatC Total 4 of points f iltered=' , 16) 

! Release classmap memory: 
Call Free<pcls) 
Call Free (pcls2) 

If (inflg .eq. 1) then . 

Call Free(pim2) 'release unneeded image memory 

call Irgb_write_isl(outf il, iml, nbx r nbands, npixm r nlinm 

npix f nlin, pixf ormo, 0) 

Else ^ ^ , 

Call Free(piral) ! released unneeded image memory 

call Irgb_write_isl (outfil, im2, nbx, nbands, npixm, nlinm 

npix, nlin, pixf ormo, 0) 

Endif 

deltime=Dtime (timearray) 
write (6, 600) timearray (1) , timearray (2 ) 
format^/' Buster Erode__Am: : user time (sec):'/ 

f8.2,'; system time (sec) : ' , f 8.2/) 

write (6, 850) outfil 

formate Cleaned output image file is ',a50) 
end 
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c Binput_isl : (10-7-92) Reads one band of an isl-format image into 
c a single byte-type array. Input is clipped to range 0-255. 



30 



35 



50 



55 



senses 



Subroutine Binput isl (inf il, im, npixm, nlinm, npix, nlin f pixf orm, 
1 band,nbanHs) 

5 

Character*50 infil 
Character cval 

Byte im (npixm, nlinm) ,bval 

Integer*2 ibuff (4096) , npix,nlin 

10 

integer imgin, ostatus, gstatus 
integer opnimg, getdef , rdline 
integer nbands, pixf orm, band, iband, ival 
integer rstatus , npx, nln, npixm, nlinm 
integer arrtyp, lininc, jtype 

'5 Equivalence (cval,bval) 

i — — 

c include ' /local/include/iopackagef ' 

!Above comment-out must be removed after f77cvt conversion!!!! 

arrtyp-IDINT2 

20 c . . . Open input image : 

imgin*0 

ostatus=opnimg (imgin, infil, 0, .false. ) 
If (ostatus .ne. 5Y5NRM) then 

write (6, 10) infil 
10 formate In image ',a50,' not open-abort') 

Call Exit(l) 

25 . Endif 

gstatus=getdef (imgin, npx, nln, nbands, pixf orm) 
If (gstatus -ne. 5YSNRM) then 
write (6,20) 

20 formate Gstatus error on input — abort') 

Call Exit(l) 

Endif 

npix^npx 
nlin^nln 
write (6, 45) 

45 format (/' "Binput isl": Npix Nlin Type Nbands Band') 

write (6, 5TJ) npix, nlin, pixf orm, nbands, band 
50 format (14x,2i5,i6, i8, i5) 

c.Read each line of input image: 
lininc=nlin/4 
jtype«lininc 
rstatus«0 

40 Do j=l,nlin 

c ... Read input lines & calc. sum: 

line-j-1 
iband-band -1 

rstatus=rdline (imgin, line, iband, ibuff, 4 096, arrtyp) 
If (rstatus .ne. SYSNRM) then 
45 write (7, 80) j 

80 formatC Error reading line -abort') 

Call Exit (1) 
Endif 



Do i*l,npix 
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End do 



Enddo 



ival-MaxO <ibuf f (i) , 0) 
ival=MinO <ival,255) 
cval=Char (ival) 
j)=bval 



!End line loop! 



.Close images 

Call clsimg (imgin) 

return 

end 
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Dedirt5x5_byte2_nopaa: 9 
Eliminates "dead zone" 

of 2 pixels around image edges (i.e. cleaning now takes place at all 
pixels of the image) . Only need one clean pixel in a 5x5 neighborhood t 
clean a dirty border pixel (i.e. nmin_edge=l) - 

Output non-dirt or cleaned classes in classmap have a value of "0". 
Hote that within one call of Dedirt, cleaned GLs are *not* used 
to clean their neighbors; this will, however, occur in 
subsequent iterations (further calls to dedirt) . 

Ad hoc "dirt* remover. Assumes that dirt is detectable as 
a simple grey-level threshold in the input classmap. 
Filters all pixels *above* a certain 
threshold grey level ("gth") . 

Note that a 1-pixel dilation of the classmap is suggested prior 
to first using this subroutine, to avoid propagation of dirt 
"ghosts" . 

Image arrays: (Max dims : npixm, nlinm; Actual dims: npix,nlin) 
iral: Input (dirty) 3-band image 
els: Input dirt 1-band classmap 
cls2: Output dirt 1-band classmap 
lm2: Output (cleaned) 3-band image 

Input: np,nl :Size of filter window. 

dgl : Threshold classmap GL value req'd for cleaning, 

fmin : Fraction of window needed to be "clean 0 to be filtered 

(l=<nmin<npMnl-l) /2 ). 

Output :cnt :No. of pixels which were filtered. 



Subroutine Dedirt 5x5 Byte2_Nopad (iml, els, im2, cls2, npixm, nlinm, npix, nlir. 
1 dgl, cnt, fminT 

Integer*2 iml (nbx , npixm, nlinm) , ixn2 (nbx, npixm, nlinm) 

Integer *2 npix,nlin 

Integer npixm, nlinm, np,nl,np2,nl2, cnt, dent, i0,ix, jO, jx,dgl 
Integer nbx, nb, num__neigh, nmin, delx, dely 
Integer nmin_edge , nmin_cent 

Byte els (npixm, nlinm) , cls2 (npixm, nlinm) ,bdgl 

Real*4 sum (3) , wtsum, wt (3, 3) , fmin 
!fmin s =.24 implies nmin=6 for 5x5 window. 

Character cdgl 

Data np,nl/5,5/ 'Erosion window si2e! 

Data wt(l,l),wt(2,l),wt(3,l)/l.,l.,0.5/ !Wts specific tr 

Data wt (1,2) , wt (2, 2) ,wt (3,2) /I., 0.707,0.4472/ 
Data wt (1,3) ,wt(2,3) ,wt (3, 3) /0 .5, 0 . 4472, 0 . 3536/ 
Data nmin_edge/0/ 

Equivalence (cdgl,bdgl) 

cnt-0 
dcnt«0 
nl2=nl/2 
np2«np/2 

nmin_cent=Nint (fmin^np^nl) 
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cdgl=Char (dgl) !0 
write (6, 10) 

formate Dedirt5x5 Byte2_Nopad: Npix Nlin Np Nl Fmin Nmin Dgl' ) 
write (6, 11) npix, nTin, np,nl, fmin, nmin cent, dgl 



10 

11 format <23x, 2i5, 2i4, f 8 . 3, i5, i6/) 



! Slide window across image: 
Do j=l,nlin 
Do iol,npix 
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!i, j are coordinates of center of window! 
!Test local value of classmap: 

If <cls(i, j) .eq. bdgl) then 'Filter the point 
dcnt^dcnt +1 
j0-j-nl2 
jx=j+nl2 
j0«Max0( jO, 1) 
1x«Min0< jx,nlin) 
i0=i-np2 
ix«i+np2 
i0=Max0U0, 1) 
ix=MinO (ix, npix) 

!iO,ix,jO,jx are coordinates of extrema of 

! current window 
Do ib=l,nb 

sum(ib) =s0. 

Enddo 

wtsum=0- 

nuro_neigh=0 

!Calc. over local "dirty* window: 
Do j2=j0,jx 

Do i2=i0,ix 

!If neighbor is *not* dirt, then: 
If (cls(i2,j2) -ne. bdgl) then 

ICalc. a distance-weight for that GL: 

!dis= (Float <i -i2))**2 + (Float (j -j2))**2 

!dis=Scsrt (dis) 

!wt*l./dis 

! wtsum=wtsum +wt 

delx=Int (Abs (i -i2) +1.) 
dely=Int (Abs(j -j2) +1.) 
wtsum=wtsum +wt (delx, dely) 

Do ib«l,nb 

sum(ib)»sum(ib) +wt (delx, dely) tt iml(ib,i. 

Enddo 

! Count no. of non-dirt window pixels: 
num neigh=num neigh +1 
Endif ~ 

Enddo 
Enddo 

JCalc. local value of "nmin*': 

If (j .le. nl2 .or. j .ge. (nlin -nl2 +1)) then 
nrain«nmin_edge !Top or bottom border! 

Else if <i .le. np2 -or. i .ge. (npix -np2 +1)) then 
nmin«nmin - edge !Left or right border! 

Else 

nmin°nmin_cent !Not in a border region! 
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ndif 11 

lit have min. no. of non-dirt neighbors, then 
! replace dirty GL w/ wted ave of nondirty GLs: 
If (num neigh .gt. nmin) then JReplace the center GL: 
5 Do , nb 

im2 (ib, i, j > =Nint (sura (ib) /wtsum) 

Enddo 

cls2(i,j)=0 iDon't filter point again! 
cnt°cnt +1 

Else !Not enough non-dirt pixels within window! 
70 !Keep the old (dirty) GL for now: 

Do ib=l,nb 

im2(ib,i, j)=iml (ib,i, j) 

Enddo 

cls2 (i, j)=cls(i, j) 
Endif 

75 Else ! Point is non-dirt, no interpolation! 

!Copy input GLs to output: 
Do ib=l,nb 

im2 (ib, i, j)~iml (ib, i, j) 

Enddo 

cls2(i, j)°cls(i, j) 

Endif 

20 Enddo !End °i° loop 

Enddo !End ° loop 

write (€,500) dcnt,cnt 
500 foraatC No. of dirt pels:',i6r'; no. of cleaned pels:',i6) 

If (dent .ot. 0 .and. cnt .eq. 0) then t 
? , write (6, «>' Warning— stable w/ some dirt remaining!!!' 

write (6, «M' [aborting]' 
Call Exit(l) 

Endif 

return 

end 
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13 

c irgb_write_isl : Writes a 3-band image from a Int*2 array. 

Subroutine Irgb write isl (bname,b, nbx, nb, npixm, nlinm, npix, nlin, 
s l pixf orm, lunprD 

Character*50 bname,outfil 

Integer* 2 b(nbx, npixm, nlinm) , npix, nlin, ibuf f (4096, 3) 
integer imgout,pstatus, ostatus 
integer nbands, pixf orm, nbx, nb, npixm, nlinm 
integer wstatus, npx, nln, ival 
integer arrtyp, lininc, jtype,pixdira, lunprt 
10 integer putdef , opnimg, wrline 

Data pixdim/4096/ 

include ' /local/include/ iopackagef ' 

arrtyp=IDINT2 
ribands-rib 

75 

write(6,45) npix, nlin, nbands , pixf orm t 
45 format <2x, ' w Irgb_write_isl" : : Npix:', 15,'; Nlin: ',15,'; Nb:',a2, 

1 '; Type:',i2) 

If (lunprt .gt. 0) then 

write (lunprt, 45) npix, nlin, nbands, pixf orm 

c...Open output image: 

imgout=0 
npx=npix 
nln=nlin 
outfil=bname 

25 pstatus=putdef (imgout, npx, nln, nbands, pixf orm) 

If (pstatus .ne. SYSNRM) then 

write (6, *)' Pstatus error for output— abort 
stop 

Endif 

ostatus=opnimg (imgout, outfil, 1, .false. ) for write-only! 

If (ostatus -ne. SYSNRM) then 
30 write(6,*)' Output image not open — abort' 

stop 

Endif 

c... Write each line of output image: 
lininc=nlin/5 
jtype=lininc 
35 Do j=l,nlin 

c ... load 2-d array & write lines: 

Do i=l,npix 

Do ib=l, nbands 

ibuff (i,ib)=b(ib,i, j) 

Enddo 
Enddo 

40 If (pixf orm .eq. IDBYTE) then 

Do ib=l, nbands 
Do i=l,npix 

ival=MaxO (ibuff (i,ib) , 0) 
ibuff £i, ib) *MinO (ival, 255) 

Enddo 
Enddo 

0 Else if (pixform .eq. IDINT1) then 

Do ib«l, nbands 
Do i=l,npix 

ival«MaxO (ibuff (i,ib) , -128) 
ibuff (i, ib) «MinO (ival, 127) 

Enddo 

50 Enddo 
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* 4 

Else if (pixrorm .eq. IDINT2) tten 
continue 

Else if (pixform .eq. IDREA4) then 
continue 

Else 

write (6, 100} pixform _ 
100 formate Pixf ormo=' , 13, ' not implemented—stop ) 

stop 

Endif 

line=j -1 , 
*rstatus=wrline (imgout , line, -1, ibuf f , paxdxm, arrtyp) 
If (wstatus -ne. SYSNRM) then 

write (6,*)' Error in wstatus — abort' 

stop 

Endif 

c If < j -ge. jtype) then 

c write<6,200) j l , x 

c 200 format UOx, ' Line ',±4,' has been written ) 

c jtype=jtype +lininc 

c Endif 

Enddo I End line loop! 

c. Close image: 

Call els img (imgout ) 

return 

end 



25 

Claims 

1. A method for the detection and removal of defects in digital images comprising the steps of: 

a) preprocessing a digital image represented by pixel values to form a feature image in which defect 
30 pixels have large local contrast from their neighboring non-defect pixels; 

b) testing the value of each feature image pixel value to determine if the value is within a range of 
pixel values expected of defects; 

c) edge-preserving spatial filtering of each of the feature image pixel values that are within the range 
of expected defects; 

35 d) forming a residual value for each pixel as a function of the difference between the corresponding 

feature image pixel value and the filtered value of step c); 

e) testing each residual value of step d) to determine if the residual value is within a range of 
. residuals values expected of defects; 

f) forming a map in which each pixel that is within the range of expected defects according to steps 
40 b) and e) is marked as a defect; and 

g) correcting the digital image using the map created in step f). 

2. A method for the detection and removal of defects in digital images comprising the steps of: 

a) preprocessing a digital image represented by pixel values to form a feature image in which defect 
45 pixels have large local contrast from their neighboring non-defect pixels; 

b) testing the value of each feature image pixel value to determine if the value is within a range of 
pixel values expected of defects; 

c) edge-preserving spatial filtering of each of the feature image pixel values that are within the range 
of expected defects; 

50 d) forming a residuals image which consists of the difference between the feature image pixel values 

and the filtered values of step c); 

e) testing each pixel value of the residuals image of step d) to determine if the residual value is 
within a range of residuals values expected of defects; 

f) forming a map in which each pixel that is within the range of expected defects according to steps 
55 b) and e) is marked as a defect; and 

g) correcting the digital image using the map created in step f), 

3. A method for the detection and removal of defects in digital images comprising the steps of: 
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a) preprocessing a digital image represented by pixel values to form a feature image in which defect 
pixels have large local contrast from their neighboring non-defect pixels; 

b) testing the value of each feature image pixel value to determine if the value is within a range of 
pixel values expected of defects; 

5 c) edge-preserving spatial filtering of each of the feature image pixel values; 

d) forming a residual value for each pixel as a function of the difference between the corresponding 
feature image pixel value and the filtered value of step c); 

e) testing each residual value of step d) to determine if the residual value is within a range of 
residuals values expected of defects; 

to f) forming a map in which each pixel that is within the range of expected defects according to steps 

b) and e) is marked as a defect; and 

g) correcting the digital image using the map created in step f). 

4. The method according to Claim 1, 2, or 3, wherein the map formed in step f) is conditioned by 
15 morphological dilation of the pixels marked as defect. 

5. The method according to Claim 1, 2, or 3, wherein the map formed in step f) is conditioned by 
morphological erosion of the pixels not marked as defect. 

20 6. The method according to Claim 1, 2, or 3, wherein the map of step f) is formed with each pixel marked 
according to the magnitude of its values according to steps b) and e). 

7. A method for the detection and removal of defects in digital images represented by multiple spatial 
resolutions ordered from lower to higher, comprising the steps of: 

25 a) preprocessing a lower spatial resolution representation of the digital image by forming a feature 

image in which defect pixels have large local contrast from their neighboring non-defect pixels; 

b) testing the value of each feature image pixel value to determine if the value is within a range of 
pixel values expected of defects; 

c) edge-preserving spatial filtering of each of the feature image pixel values that are within the range 
30 of expected defects; 

d) forming a residual value for each pixel as a function of the difference between the corresponding 
feature image pixel value and the filtered value of step c); 

e) testing each residual value of step d) to determine if the residual value is within a range of 
residuals values expected of defects; 

35 f) forming a map in which each pixel that is within the range of expected defects according to steps 

b) and e) is marked as a defect; 

g) correcting the lower resolution representation of the digital image using the map created in step 

f) ; and 

h) correcting the higher resolution representations of the digital image using the map created in step 

40 f). 

8. An apparatus for the detection and removal of defects in digital images comprising: 

means for preprocessing a digital image represented by pixel values to form a feature image in 
which defect pixels have large local contrast from their neighboring non-defect pixels; 
45 first means for testing the value of each feature image pixel value to determine if the value is within 

a range of pixel values expected of defects; 

filter means for edge-preserving spatial filtering of each of the feature image pixel values; 

means for forming a residual value for each pixel as a function of the difference between the 
corresponding feature image pixel value and the filtered value from said filter means; 
so second means for testing each residual value from said forming means to determine if the residual 

value is within a range of residuals values expected of defects; 

means for forming a map in which each pixel that is within the range of expected defects according 
to said first and said second testing means is marked as a defect; and 

means for correcting the digital image as a function of the formed map. 

55 

9. An apparatus for the detection and removal of defects in digital images represented by multiple spatial 
resolutions ordered from lower to higher, comprising: 

means for preprocessing a lower spatial resolution representation of the digital image by forming a 
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feature image in which defect pixels have large local contrast from their neighboring non-defect pixels; 

first testing means for testing the value of each feature image pixel value to determine if the value 
is within a range of pixel values expected of defects; 

filter means for edge-preserving spatial filtering of each of the feature image pixel values that are 
s within the range of expected defects; 

residual means for forming a residual value for each pixel as a function of the difference between 
the corresponding feature image pixel value and the filtered value of said filter means; 

second testing means for testing each residual value from said residual means to determine if the 
residual value is within a range of residuals values expected of defects; 
10 mapping means for forming a map in which each pixel that is within the range of expected defects 

from said first and said second testing means is marked as a defect; 

means for correcting the lower resolution representation of the digital image using the map created 
by said mapping means; and 

means for correcting the higher resolution representations of the digital image using the map 
75 created by said mapping means. 

10. A method for the detection and removal of local defects in digital images comprising the steps of: 

a) generating a feature image from an original digital image; 

b) creating an EPS-residuals image using the feature image; 

20 c) testing and mapping image pixels based upon the values of the EPS residuals image and the 

pixel values of the feature image; 

d) modifying the mapped image as a function of mapped values of neighboring pixels; 
■ e) cleaning of defect pixels by an EPS-substitution process or by a greyscale erosion process; and 
f) cleaning of defect pixels in different spatial resolution representations of the original image if 
25 different spatial resolution representations are present. 

11. A method for the detection and removal of defects in digital images comprising the steps of: 

a) preprocessing a digital image to form a feature image; and 

b) testing each feature image to determine if a value is within a range of pixel values expected of 
30 defects; 

c) edge-preserving spatial filtering of each of the feature image pixel values that are within the range 
of expected defects. 
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